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Executive  Summary 


This  final  report  presents  the  results  of  the  research  performed  under  ONR  grant  number  N0001 4-07-1- 
0043  over  the  period  of  February  1st,  2007  to  January  1st,  2010.  The  research  team  working  on  this 
project  consists  of  Prof.  Moeness  Amin  (PI),  Prof  Fauzia  Ahmad  (Research  Professor),  Dr.  Yeo-Sun 
Yoon  (Postdoctoral  Fellow),  and  Mr.  Pawan  Setlur  (Graduate  Student).  The  research  efforts  over  the  life 
of  this  grant  have  evolved  around  1)  Maximum  Likelihood  and  Suboptimal  Schemes  for  Micro-Doppler 
Estimation  using  Carrier  Diverse  Doppler  Radars,  2)  Dual  Frequency  Doppler  Radars  for  Indoor  Range 
Estimation:  Cramer-Rao  Bound  Analysis,  3)  Optimal  Waveform  Design  for  Improved  Indoor  Target 
Detection  in  Sensing  Through-the-Wall  Applications,  4)  Matched-Illumination  Waveform  Design  for  a 
Multistatic  Through-the-Wall  Radar  System,  5)  BS-MUSIC  for  High  Resolution  Imaging  in  Through-the- 
Wall  Radar  Imaging  Applications. 

Each  one  of  the  above  contributions  forms  a  chapter  of  this  report.  Each  chapter  has  its  own  Ab¬ 
stract,  Introduction,  Conclusion,  and  References  It  also  has  its  own  equation  and  figure  numbers. 

1.  Maximum  Likelihood  and  Suboptimal  Schemes  for  Micro-Doppler  Estima¬ 

tion  using  Carrier  Diverse  Doppler  Radars 

Carrier  diverse  radars,  known  as  dual  frequency  radars,  employ  two  different  frequencies,  and  can  be 
effective  in  determining  the  moving  target  range  in  urban  sensing  and  through-the-wall  radar  applications. 
In  this  report,  we  derive  the  maximum  likelihood  (ML)  estimator  for  the  micro-Doppler  motion  parame¬ 
ters  from  the  dual  frequency  radar  returns.  Micro-Doppler  signatures,  which  are  commonly  associated 
with  vibrating,  oscillating,  and  rotating  objects,  have  emerged  to  be  an  important  tool  in  target  detection 
and  classification.  Unlike  linear  models,  the  respective  ML  estimator  does  not  assume  a  closed  form 
expression  We  solve  the  ML  estimator  for  dual-frequency  radar  operations  by  applying  an  iteratively 
reweighted  nonlinear  least  squares  algorithm  (IRNLS),  which  is  initiated  using  suboptimal  solutions.  The 
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ML-IRNLS  algorithm  is  applied  to  both  simulated  and  experimental  radar  returns  for  estimating  the  range 
and  the  motion  parameters  of  indoor  targets. 

2.  Dual  Frequency  Doppler  Radars  for  Indoor  Range  Estimation:  Cramer- 

Rao  Bound  Analysis 

Single  frequency  Doppler  radars  are  effective  in  distinguishing  moving  targets  from  stationary  targets,  but 
suffer  from  inherent  range  ambiguity.  With  a  dual-frequency  operation,  a  second  carrier  frequency  is 
utilized  to  overcome  the  range  ambiguity  problem.  In  urban  sensing  applications,  the  dual-frequency 
approach  offers  the  benefit  of  reduced  complexity,  fast  computation  time,  and  real  time  target  tracking. 
We  consider  a  single  moving  target  with  three  commonly  exhibited  indoor  motion  profiles,  namely, 
constant  velocity  motion,  accelerating  target  motion,  and  micro-Doppler  motion.  RF  signatures  of  indoor 
inanimate  objects,  such  as  fans,  vibrating  machineries,  and  clock  pendulums,  are  characterized  by  micro- 
Doppler  motion,  whereas  animate  translation  movements  produce  linear  FM  Doppler.  In  this  report,  we 
derive  Cramer-Rao  bounds  (CRB)  for  the  parameters  defining  indoor  target  motions  under  dual-frequency 
implementations.  Experimental  data  is  used  to  estimate  micro-Doppler  parameters  and  to  validate  the 
CRBs. 

3.  Optimal  Waveform  Design  for  Improved  Indoor  Target  Detection  in  Sens¬ 

ing  Through-the-Wall  Applications 

The  report  also  deals  with  waveform  design  for  improved  detection  and  classification  of  targets  behind 
walls  and  enclosed  structures.  The  target  impulse  response  is  incorporated  in  an  optimum  design  of  the 
transmitted  waveform  which  aims  at  maximizing  the  signal-to-interference  and  noise  ratio  (SINR)  at  the 
receiver  output.  The  interference  represents  signal-dependent  clutter  which,  along  with  the  wall,  degrades 
the  receiver  performance  compared  to  the  free-space  and  zero-clutter  case.  Computer  simulations  show 
sensitivity  of  the  optimum  waveform  to  target  orientation  but  depict  an  SINR  enhancement  over  chirped 
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waveform  radar  emissions  at  all  aspect  angles.  Numerical  electromagnetic  modeling  is  used  to  provide  the 
impulse  response  of  typical  indoor  stationary  targets,  namely,  tables,  chairs,  and  humans. 


4.  Matched-Illumination  Waveform  Design  for  a  Multistatic  Through-the- 

Wall  Radar  System 

We  present  the  matched  illumination  waveform  design  for  improved  target  detection  in  through-the-wall 
radar  imaging  and  sensing  applications.  We  consider  a  multistatic  radar  system  for  detection  of  stationary 
targets  with  known  impulse  responses  behind  walls.  The  stationary  and  slowly  moving  nature  of  typical 
indoor  targets  relaxes  the  orthogonality  requirement  on  the  waveforms,  thereby  allowing  sequential 
transmissions  from  each  transmitter  with  simultaneous  reception  at  multiple  receivers.  The  generalization 
of  the  matched  illumination  waveform  design  concept  from  a  monostatic  to  a  multistatic  setting  casts  the 
indoor  radar  sensing  problem  in  terms  of  multiple-input  multiple-output  (MIMO)  operations  and  puts  in 
context  the  offering  of  MIMO  to  urban  sensing  and  imaging  of  targets  in  enclosed  structures.  Numerical 
electromagnetic  modeling  is  used  to  provide  the  impulse  response  of  typical  behind-the-wall  stationary 
targets,  namely  tables  and  humans,  for  different  target  orientations  and  at  various  incident  and  reflection 
angles.  Simulation  results  depict  an  improvement  in  the  signal-to-clutter-and-noise-ratio  (SCNR)  at  the 
output  of  the  matched  filter  receiver  for  multistatic  radar  as  compared  to  monostatic  operation. 

5.  BS-MUSIC  for  High  Resolution  Imaging  in  Through-the-Wall  Radar  Imag¬ 

ing  Applications. 

The  MUSIC  algorithm  is  a  high-resolution  direction  finding  technique  which  has  been  successfully 
applied  to  enhance  radar  imaging  in  inverse  synthetic  aperture  radar  (ISAR).  Although  this  signal  sub¬ 
space-based  method  has  proven  effective  when  dealing  with  point  targets  and  high  SNR,  it  may  fail  to 
work  when  directly  applied  to  extended  targets  or  target  returns  of  low  SNR.  The  Beamspace  MUSIC 
(BS-MUSIC),  in  which  the  MUSIC  algorithm  is  applied  to  multiple  beams  is  capable  of  locating  spatially 
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extended  targets  in  low  SNR  environments.  We  propose  BS-MUSIC  as  an  image  formation  method  for 
indoor  radar  imaging  problems  and  sensing  through-the-wall  applications.  We  compare  BS-MUSIC 
performance  to  conventional  beamforming  and  element-space  MUSIC.  Imaging  results  with  both  synthe¬ 
sized  and  real  data  demonstrate  the  advantages  of  the  proposed  algorithm  in  depicting  targets  behind 
walls. 
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1.  Maximum  Likelihood  and  Suboptimal  Schemes  for  Micro-Doppler 
Estimation  using  Carrier  Diverse  Doppler  Radars 


Abstract 

Carrier  diverse  radars,  known  as  dual  frequency  radars,  employ  two  different  frequencies,  and  can  be 
effective  in  determining  the  moving  target  range  in  urban  sensing  and  through-the-wall  radar  applications. 
In  this  chapter,  we  derive  the  maximum  likelihood  (ML)  estimator  for  the  micro-Doppler  motion  parame¬ 
ters  from  the  dual  frequency  radar  returns.  Micro-Doppler  signatures,  which  are  commonly  associated 
with  vibrating,  oscillating,  and  rotating  objects,  have  emerged  to  be  an  important  tool  in  target  detection 
and  classification.  Unlike  linear  models,  the  respective  ML  estimator  does  not  assume  a  closed  form 
expression.  We  solve  the  ML  estimator  for  dual-frequency  radar  operations  by  applying  an  iteratively 
reweighted  nonlinear  least  squares  algorithm  (IRNLS),  which  is  initiated  using  suboptimal  solutions.  The 
ML-IRNLS  algorithm  is  applied  to  both  simulated  and  experimental  radar  returns  for  estimating  the  range 
and  the  motion  parameters  of  indoor  targets. 
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1.1  Introduction 


Urban  sensing  and  through-the-wall  radar  imaging  address  the  desire  to  detect,  locate,  and  classify  both 
animate  and  inanimate  targets  [l]-[4].  Range  estimation  is  typically  performed  by  linear  frequency 
modulated  radars,  pulse  compression  radars,  or  pulse  Doppler  radars  [5],  Such  radar  systems  are  wide¬ 
band  so  as  to  meet  the  range  resolution  requirements.  However,  the  operational  logistics  and  system 
requirements  for  urban  sensing,  such  as  cost,  hardware  complexity,  and  portability,  may  impede  or 
prohibit  the  use  of  such  radar  systems.  Further,  bandwidth  allocation  issues  may  arise  since  radio  fre¬ 
quency  (RF)  penetration  through  the  walls  follows  a  lowpass  filtering  model  with  typical  cutoff  in  the  low 
GHz  range,  where  much  of  the  RF  spectrum  may  be  jammed  or  taken  over,  in  part  or  fully  by  adversaries 
or  other  emitters.  On  the  other  hand,  a  dual-frequency  approach  for  target  range  estimation,  combined 
with  wide  array  aperture,  can  meet  the  requirements  of  different  system  operation  modes,  and  is  likely  to 
emerge  as  the  preferred  approach  in  most  urban  sensing  and  rescue  missions  [1],[6].  The  dual  frequency 
or  carrier  diversity  is  induced  by  using  two  different  carrier  frequencies  which  are  selected  to  achieve  a 
desirable  maximum  unambiguous  range.  The  latter  is  important  to  allow  a  unique  range  estimate  of  a 
target,  and  should  be  based  on  the  a  priori  knowledge  (possible  through  aerial  mapping  or  ground  access) 
of  the  spatial  extent  of  the  urban  structure  under  surveillance.  The  technique  of  employing  two  frequen¬ 
cies  to  estimate  range  has  been  used  in  many  other  radar  applications  [7],  [8]. 

In  this  chapter,  we  apply  a  dual  frequency  radar  for  target  range  and  parameter  estimation.  We  con¬ 
sider  the  micro-Doppler  (MD)  target  motion  profile  [9].  Micro-Doppler  analysis  has  been  used  in  many 
applications  for  human  gait  analysis,  multistatic  radar  applications,  etc.,  such  works  can  be  seen  for 
example  in  [10],  [11],  and  references  therein.  RF  signatures  of  indoor  inanimate  objects,  such  as  fans, 
vibrating  machineries,  and  clock  pendulums,  and  animate  objects,  like  the  limbs  in  human  gait  are 
characterized  by  MD  motion.  Translational  motions,  producing  constant  velocity  or  accelerating  velocity, 
respectively  produce  complex  sinusoids  and  linear  frequency  modulation  to  the  incident  waveform.  The 
ML  techniques  for  parameter  estimation  of  such  returns  has  been  treated  in  [12],  [13],  and  references 
therein.  However,  the  ML  for  micro  Doppler,  which  gives  rise  to  sinusoidal  FM  signals,  has  not  yet  been 
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examined  and  is  the  subject  of  this  chapter.  We  derive  the  maximum  likelihood  optimal  estimator  for  MD 
motion  parameters.  The  results  are  compared  with  the  Cramer-Rao  lower  bounds  (CRBs)  for  MD  motion 
parameters  which  were  recently  derived  in  [14]  for  dual-frequency  operations. 

We  consider  a  single  moving  target  whose  MD  motion  profile  can  be  modeled  by  a  finite  number  of 
parameters.  The  MD  is  further  classified  as  rotational  or  vibrational  MD  based  on  radar  cross  section 
(RCS)  fluctuations.  Maximum  likelihood  (ML)  technique  for  motion  parameter  estimation  is  then 
formulated  and  solved  using  step  wise  concentration  to  obtain  an  iteratively  reweighted  least  squares 
algorithm.  The  iterative  algorithm  is  initialized  using  suboptimal  estimates  and  applied  to  real  radar 
returns  to  obtain  ML  estimates  of  the  MD  parameters.  It  is  noted  that  the  focus  of  this  chapter  is  on  single 
antenna  operation.  For  a  multiple  target  scene,  assuming  that  the  targets  are  separable  in  cross-range  and 
spatial  processing  (beamforming)  is  used  in  conjunction  with  the  dual-frequency  radar,  the  ML  analysis 
presented  in  this  chapter  is  applicable  to  each  individual  separated  target  return. 

A  brief  outline  of  the  chapter  is  as  follows.  Section  1.2  describes  the  signal  model.  In  Sections  1.3 
and  1.4  we  discuss,  respectively,  the  ML  and  suboptimal  estimation  schemes  for  the  micro-Doppler 
motions.  Section  1.5  contains  the  simulation  and  experimental  results,  followed  by  the  conclusions  in 
Section  1.6. 


1.2  Signal  Model 

The  signal  returns  for  the  dual  frequency  Doppler  radar  after  down  conversion  to  baseband  and  using  N 
samples  are  given  by, 


jc1(n)  =  51(/i)  +  vi(n)  =  Ai(/i;\j»')exp|  j—  —  R(n,y)  +vi(n),  /  =  !,  2,  and n  =  0,l,...,Af-l 


V  c  j  (1) 

Ejv,(n)v2*(£)J  =  0  Vn,£;  £|vJ(n)vI*(A:)j  =  0  i  =  \,2,E{vi(n)v‘(k'^  =  <jf  \/n  =  k,i  =  1,2 

where  f:,  i  =  1,2  are  the  carrier  frequencies,  c  is  the  speed  of  light  in  free  space,  and  the  superscript  ‘  *  ’ 
denotes  the  complex  conjugate  operation.  The  target  range,  /?(/j;\|/),  is  parameterized  by  a  vector 

vpe  of  P  desired  parameters.  The  amplitude  («;»}/') ,  measured  at  the  t-th  frequency,  captures  the 


17 


time  variations  in  the  target  RCS,  and  is  a  function  of  a  subset  \\f'  of  the  parameter  vector  y.  The  noise, 
VjO  and  v>20,  at  the  two  carrier  frequencies,  are  assumed  to  be  complex  circular  AWGN  and  uncorre¬ 
lated.  Further,  the  noise  sequences  are  i.i.d  for  each  carrier  frequency.  Note  that  the  model  in  (1) 
corresponds  to  the  discrete-time  equivalent  of  the  continuous-time  signals,  i.e.,  x :,.(«)  =  Jtf(/f  7^),  where 

Ts  is  the  sampling  period  and  has  been  suppressed  in  the  notations  for  convenience.  The  returns  in  ( 1 )  can 

be  statistically  characterized  by  a  multivariate  complex  Gaussian  probability  density  function  (pdf), 
p(x;s). 


p(x\s)  =  - 


1 


exp(-(x-p)wC  '(x-ji)) 


n2N  det(C) 


where  the  received  signals  at  the  two  frequencies  are  appended  to  form  a  long  vector, 

X  =  s  +  v  =  [x1(0),X|(1),...,xi(jV-1),x2(0),x2(1),...,x2(jV  -l)f  =[x,  x2]T 
s  =  [s](0),s](l),...,sl(N -\),s2(0),s2(\),...,s2(N -l)]T  =[s}  s  J 
v  =  [ v,  (0),  v,  (1), ....  v,  (N  - 1),  v2  (0),  v2  (1), . . . ,  v2  ( A/  -  l)]r  =  [  v,  v2f 


(2) 


(3) 


The  mean  \i  =  £{x} equals  [s}  s2]7  and  the  covariance  matrix  Cis  Hermitian  with  the  following  diagon¬ 

al  structure, 


„  crrl  0VvA/1 
C  =  £{(x-p)(x-p)w}=  1  NxN 

®NxN  J 


(4) 


where  I  is  an  identity  matrix  of  dimensions  NxN. 


1.3  Maximum  Likelihood  Estimation 

We  consider  the  noise  free  return  s,  which  is  a  function  of  R(n,v\>),  /r((n;i)/),  and  ft,  (  =  1,2,  to  be 

parameterized  by  the  vector  of  P  desired  parameters.  Hereafter,  for  notational  succinctness,  we  denote 

s  =  s(\p) .  The  ML  estimator  for  \j/  =  [yfx,  \f/2,  y/2 . i//P]T  is  defined  as  [15], 

»p  =  argmax/r(x;5(v}/))  =  argmax  ln(/r(x;.s(v|/)))  (5) 

v  v 
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If  the  elements  of  the  covariance  matrix,  i.e.,  <72,  /  =  1,  2,  are  unknown,  then  using  (2)  and  ignoring  the 
constant  terms,  the  ML  estimator  of  the  complete  parameter  vector,  0  =  [\|/7  ,rr2,  <j\]‘  can  be  expressed 


as, 


( 


0  =  arg  max 
o 


-N  ln((X2 )  -  N  In (oi )  - J 


-S2(V|/)|| 


2  "N 


(6) 


Computing  the  log  likelihood  (LL)  score  with  respect  to  the  noise  nuisance  parameters  (J2,  i  =  1,  2,  and 
equating  it  to  zero,  we  obtain 

,2  _  IK  -Mv)|| 


N 


,  for  i  =  l,  2 


(7) 


Concentrating  the  LL  with  respect  to  of  ,  i.e.,  substituting  (7)  in  (6),  yields 


\\f  =  arg  max 


—N  In 


x,-s,(v)i  jh-MHr 

N  N 


&2  =  Ik-  -s,-(N»|| 

'  N 


=  arg  min  in 


IT 


N 


2  \ 


(8a) 


(8b) 


The  ML  estimator  for  the  noise  nuisance  parameters  is  provided  in  (8b),  and  is  the  same  as  the  expression 
in  (7)  with  \]f  replaced  by  \j>.  Notice  that  the  noise  estimate  at  one  frequency  only  depends  on  the  data 
measured  at  the  same  frequency.  However,  it  also  depends  on  the  ML  parameter  estimates  obtained  from 
the  combined  frequency  information,  and  as  such,  the  problem  cannot  be  decoupled  into  two  separate  ML 
parts,  each  corresponding  to  a  single  frequency.  Equation  8(a,  b)  constitutes  the  ML  estimator  for  the  dual 
frequency  radar,  and  due  to  the  involvement  of  the  two  terms,  it  does  not  have  a  closed  form  solution.  If 
the  elements  of  the  covariance  matrix  C  are  known,  the  ML  estimator  for  \\r  takes  the  form 


(\ 


\j)  =  arg  min  (x  —  s(\j#))/7  C  J(x  -s(\|/))  =  arg  min 


X, -S|(V|/)1  |x2-s2(v)|| 


|2  \ 


'1 


(9) 
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The  step  wise  concentration  approach  can  be  applied  using  (9)  to  provide  an  iterative  solution  to 
8(a,b).  The  iterative  ML  algorithm  is  formulated  as  follows. 

1 .  Initialize  with  C0  =  I ,  where  I  is  an  identity  matrix  of  the  same  dimension  as  C0. 

2.  Using  (9),  find  the  estimate  vj/0  of  \j /. 

3.  Use  \j )0  in  (8b)  to  obtain  the  noise  variance  estimates  =||x/  —  s(-(\jjr0)|  /A,  and  construct  the 

covariance  matrix  C1  =  Diag  [a,2,  xl,^  xl],  where  the  operator  Diag[  ]  transforms  a  vector 
into  a  diagonal  matrix  and  1  =  [1,1,1..  l]lx/v. 

4.  Recursively  solve  at  the  Ath  iteration  (A  >1), 

vj )k  =  argmin  (x-  s(\\f))If  Ckl  (\-  s(\\f))y  initialized  with 
v 

^*+.)=IK-s,(V*)|2/^  Vi  =  1,2  ;  C(i+I)  :=  Diag[a?<k+l)  xl,&l<k+l)xl]  (10) 

5.  Stop  at  convergence,  or  when  an  appropriate  terminating  criterion  is  satisfied. 

The  solution  of  step  4  depends  on  the  underlying  motion  model,  as  delineated  next.  We  note  that  the 
estimate  vj>()  is  the  non-linear  least  squares  (NLS)  estimate.  Further,  in  step  4,  the  covariance  matrix  C*  is 
not  stressed  to  be  an  estimate  for  reasons  to  follow  shortly.  As  evident  from  the  above  iterative  ML 
implementation,  the  step-wise  concentration  approach  does  not  treat  the  noise  variances  as  of  known  or 
estimated  values,  but  rather  alternates  between  the  two  assumed  hypotheses.  In  other  words,  for  every 
iteration,  a  quasi-ML  objective  is  optimized.  The  step-wise  algorithm  has  been  used  in  generalized  linear 
models  in  statistical  literature  [16]  and  in  robust  statistics  for  M-estimation  [17].  It  is  often  described  as 
the  iteratively  reweighted  least  squares  (IRLS).  The  “reweighting”  at  each  iteration  occurs  in  the  estima¬ 
tion  of  \\tk  through  a  known  C*,  as  seen  in  (10).  However,  for  the  problem  at  hand,  as  shown  in  the 
following  section,  s(\| /)  is  nonlinear,  and  hence  a  more  appropriate  name  for  the  algorithm  would  be  the 

iteratively  reweighted  nonlinear  least  squares  (IRNLS).  Henceforth,  we  will  refer  to  the  iterative  ML 
algorithm  as  IRNLS.  Proof  of  convergence  of  the  IRNLS  estimates  is  provided  in  the  Appendix. 
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Three  terminating  criteria  may  be  considered  for  the  IRNLS  algorithm.  Denoting  the  negative  LL 
cost  by  an  appropriate  criterion  is  when  P(\\tk)-  where  £  is  an  arbitrarily  small  user- 

defined  positive  constant,  in  which  case  vj )k  will  be  the  final  ML  estimate.  Another  terminating  criterion 
could  be  | \\?k  “\j rk  ,|<£,  where  £  is  the  available  machine  precision  or  can  be  user-defined.  The  third 

terminating  criterion  is  simply  constraining  the  maximum  number  of  iterations,  i.e.,  iterate  (10)  until 
k  =  £max  •  We  use  this  criterion  in  the  simulations.  A  general  flowchart  of  the  ML  algorithm  is  shown  in 
Fig.  1. 

In  the  following,  we  associate  s(vj /)  to  different  MD  motion  profiles,  and  use  the  IRNLS  to  derive 
the  respective  optimal  ML  solution. 

1.3.1  Micro-Doppler 

The  MD  returns  can  be  classified  as  a)  vibrational  MD  and  b)  rotational  MD.  Although  the  phase  of  the 
returns  is  identical  in  both  cases,  a  difference  between  the  two  exists  in  terms  of  the  amplitude  or  RCS 
fluctuations. 

1.3.1. 1  Vibrational  MD 

The  vibrational  MD  arises  due  to  vibrations  of  the  scatterers  on  the  target  or  of  the  target  itself,  example 
being  a  target  moving  back  and  forth  or  undergoing  a  simple  harmonic  motion  (SHM).  The  vibrational 
MD  is  characterized  by  a  sinusoidal  instantaneous  frequency.  It  is  parameterized  by  vj/  =  [Ro,d,coo,(pa], 
where  R()  is  the  range  of  the  target  at  rest  or  simply  the  mean  range,  d  is  the  maximum  radial  displace¬ 
ment,  (Oa  is  the  vibrational  frequency,  and  <po  is  a  constant  phase.  The  time-varying  range  profile  for  this 
motion  is  given  by, 

R(n-,y)  =  Rtl  +  d  cos(a>0n  -  <po)  (11) 

Typical  indoor  vibrating  targets  have  small  displacements  relative  to  the  target  range,  especially  for 
longer  radar  standoff  distances  from  the  wall  that  are  normally  used  for  through-the-wall  radar  operations. 
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Accordingly,  the  target  aspect  angle  and  hence  the  RCS  is  considered  constant,  thereby  yielding  a 
constant  amplitude,  i.e.,  /z/(n;\|/,)  =  /7l,  /  =  1,2.  Substituting  the  expressions  for  /if. (/i; \ /)  and  R(n\y \t)  in 
(1),  we  obtain  the  signal  returns  as. 


/ 


4 xfj ( R,  +  d  cos {a)0n  -  (p„ ))  ^ 


-*;(«)  =  A  ex p  j 


+  v7(/i),  Vi  =  1,2 


02) 


v 


c 


The  radar  returns  in  (12)  represent  sinusoidal  FM  signals.  Using  the  data  vector  x,  comprised  of  the 
returns  of  (12),  in  (10)  defines  the  IRNLS-ML  estimator  for  the  vibrational  MD  motion  profile.  The 
noise-free  signal  vector  s  can  be  decomposed  into  two  multiplicative  terms,  one  containing  Ro  and  the 
other  being  a  function  of  the  remaining  parameters  in  \\i.  Accordingly,  we  express  the  IRNLS-iterations 
as 


\\tk  =argmin(x- A(v|/,)b)//CJt1(x-  A(v|//)b),  k  >  1,  initialized  with 


ai(V,):=D»expO"4tf/j.z(l)/ c),...,exp(j4nfiz(N -\)/ c)f ,  1/  :=[d,eon,<pJ 
z(n):=dcos(a)0n-<p0),  bt  ~  ptexp(j4nfiR0 /c),V/'  =  l,2 


(13) 


In  (13),  bik  denotes  the  estimate  ht  at  the  A:-th  iteration.  We  proceed  by  minimizing  (13)  with  respect  to 


b,  in  which  case  one  obtains  the  well  known  weighted  least  squares  solution,  substituting  this  into  (13) 
we  obtain  the  final  estimates  of  the  parameters  as. 
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(dk,ai,^ok)  =  argmax 

d,(o0,<p0 


N- 1 


Y  *1  (n)exp(-j4xftd  cos(  (00n-<po)l  c) 


n= 0 


2 


c Xu 


/V-l 

X  *2 («) exp(-y'4^/2 d  cos(0)on  -%)!  c) 

n=0 


'2k 


,  /:  >  1 


A*  ^  A(vi),b*  -b(y;)  =  (A»qlAi)-,A«C*lx,  ^  =^i^klc 


(14) 


A*  - 


b-. 


.  ^ 


(*+D 


x,-a,(vK)4k 


/  /V,  i  =  l,2 


In  (14),  A*  and  b^,  respectively,  denote  the  estimates  of  A  and  b  in  the  A>th  iteration,  defined  in  (13). 
Similar  notational  convention  follows  for  the  other  parameters.  The  function  maximizations  in  (14)  are 
solved  numerically.  It  is  noted  in  (14)  that  the  range  estimate  Rok  is  obtained  by  subtracting  the  mean 

phase  estimates  at  the  two  carriers  and  dividing  it  by  the  difference  of  the  carrier  frequencies.  This 
estimate  corresponds  to  the  standard  dual-frequency  range  estimate  to  avoid  the  many  ambiguous  range 
solutions  see  [7,  pg.  140].  In  other  words,  we  are  assuming  Ro  e  [0,c/  2(/2  -/,)]. 


1.3.1.2  Rotational  MD 

As  the  name  suggests,  targets  which  are  rotating  with  respect  to  a  fixed  location  radar  follow  this  model. 
Unlike  the  vibration  model,  the  signal  returns  due  to  rotation  have  RCS  fluctuations.  This  is  because  the 
radar  observes  different  elevation  aspects  of  the  target,  thereby  inducing  a  cyclic  amplitude  modulation  in 
the  radar  returns.  In  general,  the  RCS  fluctuations  are  geometry  dependent.  For  example,  the  RCS 
fluctuations  are  non-existent  for  a  rotating  sphere  since  the  sphere  is  aspect  independent,  whereas  for 
other  complex  targets,  the  RCS  fluctuates  cyclically.  In  this  chapter,  we  consider  the  sine  type  fluctuation, 
which  corresponds  well  to  a  rotating  fan  blade  [18-20].  It  must  be  noted  that  most  of  the  typical  indoor 
rotating  targets  are  fans,  which  may  either  be  ceiling  mounted,  pedestal  or  table-top.  This  is  primarily  the 
reason  for  using  the  sine  model  in  the  underlying  application  area.  As  such,  the  baseband  returns  at  the 
two  carrier  frequencies,  for  a  single  blade,  can  be  readily  shown  to  be 
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xi(n)  =  hj(rv,\ |/')exp 


•  IxJ,  Cgo  +  d  cosjcoji  -  g, )) 


(»;'/)  :=  A)'(«;i|' ')=/?,  sine 


c  y 

A 


+  v,(/i),  i  =  1,2 


4;r/,<fcos(<y„/i-pr,) 


(15) 


y 


,  v|/':=  [d,<y(,,<pj7 


The  sine  function  in  (15)  is  defined  as  sinc(jc):=sin(;c)/x  From  (15),  it  is  noted  that  the  sine  function 
has  a  cyclic  behavior  which  depends  on  the  rotational  frequency.  The  rotational  frequency,  in  turn,  is  also 
a  function  of  the  sampling  frequency,  and  hence  the  data  in  (15)  must  span  at  least  one  cycle  for  any  type 
of  processing  to  be  successfully  applied  to  it.  Extending  the  model  in  (15)  to  Q  blades,  the  return  at  the  i- 
th  frequency  is  given  by, 


x,(n)  =  YJhll)(n\\\i')exp 


q=\ 


. 4^/,  ( K  +  d  cos  -g>q)) 


+  vf.(rc),  /  =  1,2 


hiq{n\Mf')  :=  ptf^nw’)  =  piq  sine 


4yr  fjd  cosjcoji  —  gj  ’ 


(16) 


The  model  in  (16)  is  general,  in  the  sense  that  it  assumes  the  Q  blades  are  different,  and  thus  piq  is 

indexed  by  q.  If  the  blades  possess  identical  geometry,  and  no  have  manufacturing  defects,  then  as  a 
special  case  puj  =  pt,  \/q.  Further,  assuming  that  the  blades  are  placed  symmetrically  around  the  main 

rotor,  we  can  express  (pq  as 


^=«,  +  to-l)2*/G,  q  =  U  (17) 

In  this  case,  the  radar  return  from  the  rotating  object  has  a  harmonic  structure  with  harmonic  frequencies 
at  mQct)(nm  =  0,±l,...,°o  [18-20].  We  must  note  that  the  RCS,  which  is  strictly  positive  and  real,  is 
determined  from  the  magnitude  squared  of  the  noise  free  returns  in  (15)  and  (16).  The  returns  in  (16)  can 
be  rewritten  in  a  compact  vector  notation  as 
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x  =  [x[,x[]r  =A(y')b 


A(V'):= 


A"(M,,)  °  ,|1b:=[bf,bnr 

0  A22(v')1]  1  2J 


*>,  =[^pfe,2 . &,0f -  K  =Piqexp(j4xfiRo/c)’  A„(M/'):=[a,,(M/,),a,.2(v|/,),...,a,0(M/')],  i  =  I,2  (18) 

»,„(»/):=?„□  g.,,  Y„  :=[^(0;M',),ri,(l;v').-.yi?(^-Uv,)]r,  q  =  1,2,. ..,(2 

/  /x  r  j4xftd «>s( 9?  )/ c*  )/c  >4;r/,c/cos(w(,(A/-l)-®  )/r ^ 

):=[<?  ^  J 

In  (18),  the  symbol  ‘C  ’  denotes  the  Hadamard  product  or  element-wise  product.  Following  the  analysis 

for  vibrational  MD  case,  it  can  be  readily  shown  that  the  IRNLS  ML  for  rotational  MD  at  the  kth  iteration 
is  given  by 

\\f[  ^argmaxx^P^y^C^x,  initialized  with  y'  , 
y' 

P*  (*/) :=  A*  ( Af  C~k}Ak)~[  Af  C*1 
b*  =  (A?  c;‘ A*  r 1  Af  c;’x  =  [b[*  ,bT2k  f 

R,k=lTZ{b;k0b2k)x- 


(19) 


~  2 

°i{k+ 1) 


xi-w;^2  iNj=\xp,q  =  biq^q 


where  1  is  column  vector  of  dimensions  Q  x  1  whose  elements  are  comprised  of  all  ones,  and  A(\|/)  and 


b  are  defined  in  (18).  In  (19),  A*  :=  A(\j/*)  is  the  estimate  of  the  matrix  Ak.  The  function  maximiza¬ 
tions  in  (19),  as  before  is  carried  out  numerically.  A  special  arises  when  the  blades  are  all  identical.  I.e., 
P-  —  Vt/,  then  the  IRNLS  ML  derived  thus  far  still  applies,  for  example  in  (19),  we  simply  use 

Pi  =  1  jb where  the  absolute  value  is  taken  element-wise.  For  a  single  blade  return,  we  can  substi¬ 
tute  q-  I  in  (18-19)  and  proceed  with  the  analysis. 

We  use  the  CRBs  to  compare  the  mean  squared  error  (MSE)  of  the  ML  estimates  for  both  vibrational 
and  rotational  MD  motions.  The  elements  of  the  Fisher  Information  Matrix  (FLM)  for  rotational  MD 
assuming  the  sine  model  are  derived  in  Section.  III.  The  expressions  can  then  be  used  to  numerically 
evaluate  the  inverted  FIM.  The  FIM  for  vibrational  MD,  derived  in  [14],  is  a  special  case  of  the  rotational 
MD  FIM. 
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1.4  Suboptimal  Estimators 

Iterative  non-linear  schemes,  ML  or  otherwise,  require  proper  initial  estimates  for  achieving  convergence. 
The  non-linear  cost  functions  of  (8a)  and  (9)  have  multiple  local  extrema.  Therefore,  in  order  to  obtain 
\\ in  step  2  of  the  IRNLS-ML  algorithm,  initial  estimates,  obtained  from  suboptimal  estimators  which 
depend  on  the  noise  free  returns,  st(\\ /),  i  =  1,2,  can  be  used.  Below,  we  discuss  the  suboptimal  estimation 
schemes  for  both  vibrational  and  rotational  MD. 

1.4.1  Vibrational  MI) 

Ignoring  the  contribution  of  noise  for  the  time  being,  we  note  that  the  Fourier  spectrum  of  the  return, 
x /  =  1,2  in  (12)  is  not  analytic,  and  consists  of  infinitely  many  harmonics  weighted  by  Bessel 
functions  of  the  first  kind  [21], 

Xi(co)  =  2npiexp(j4nfiR0  /  c)^  jmJm(4nfid/c)5(a)-ma)0)exp(-jm<p0),  Vi'  =  1,2  (20) 

m=~oo 

where  Jm()  is  the  mh  order  Bessel  function  of  the  first  kind.  Since  the  Bessel  functions  rapidly  decrease 
in  magnitude  for  increasing  m,  the  Fourier  transform  (20)  of  the  returns  in  (12)  has  at  the  most  m  =  K 
significant  harmonics.  The  value  of  K  is  readily  inferred  from  the  spectrum.  With  this  in  mind,  we 
describe  below  the  suboptimal  estimation  procedure.  Since  the  noise  variances  are  neither  required  for 
estimating  the  parameters  suboptimally  nor  are  they  essential  to  start  the  IRLS  iterations,  their  suboptimal 
estimates  are  omitted. 

To  obtain  initial  estimates  of  0)o ,  we  choose  2AT+1  peaks  of  Xt{(D ),  i  =  l,2  and  form  the  vector 
y-  =  [Xi(O)_K)9Xi(C0LK+l)9...9Xi(0)J...9Xi(6)K_l)9Xi(6)K)Y-  The  frequencies  corresponding  to  these  peaks 
are  stored  in  a  vector  co(.  Since  the  noise  is  different  for  each  carrier  frequency,  different  perturbations  in 
the  peak  frequency  locations  can  occur.  Hence,  in  general,  o1  ^to2.  The  DC  value  in  the  Fourier  trans¬ 
form  is  important,  since  it  serves  a  reference  for  automatic  peak-picking,  and  therefore  it  must  be  included 
in  the  vector  yf..  The  suboptimal  estimates  O)oi  suhopi ,  i  =  1,2,  for  COo  can  be  obtained  as. 
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(21) 


a) 


'oi, subopt 


K' to 
(K7K) 


S  V/  =  1,  2,  K  \=[-K,-K  + 1, _ ,0, _ _  AT]Z 


Equation  (21)  is  similar  to  the  least  squares  estimator  proposed  in  [22,  eq.  (35)],  with  subtle  differences 
due  to  the  model  choice. 

In  order  to  estimate  the  remaining  parameters,  albeit  suboptimally,  one  first  needs  to  estimate  the 
amplitude  pt .  For  this  purpose,  we  employ  a  computationally  efficient  technique  using  higher  order 
statistics  which  was  proposed  in  [22,  eq.  (50-52)], 


Pi,  subopt 


=  W2 


f  j  N- 1 


\2 


N  1 


Ti ZMn)l  • v/=1’2 


(22) 


n-0 


The  above  equation  relies  on  the  higher  order  moments  of  the  Gaussian  random  variable,  and  may  not  be 
valid  for  other  types  of  distributions.  The  estimation  of  dt  depends  on  classifying  the  signal  returns  as 
wideband  or  narrowband.  We  consider  the  wideband  case1.  The  suboptimal  estimator  proposed  in  [23], 
and  also  used  in,  [22,  eq.  (38),  eq.  (40)],  and  [24]  can  be  applied  to  obtain  the  estimate  of  dt,  i  =  1,2 ,  as 


V  a  (y,  y,) 


(23) 


The  estimate  for  Ro  is  given  by, 


subopt  =  (ZX2(0)X;(0))c/47r(f2-f{)) 


(24) 


Indeed,  other  techniques,  based  on  (20),  can  be  considered  for  estimating  R().  However,  the  suboptimal 
estimates  of  (24)  were  found  to  perform  reasonably  well  in  terms  of  the  mean  square  error.  For  <po,  we 
use  the  LS  estimator  proposed  in  [22]  after  appropriately  demodulating  the  phase  estimates  for  Ro. 


^Poi, subopt 


-K7  arg(y,exp(-;4^-/.^  wfav„  /c)) 

IN2 


(25) 


1  Narrowband  case  is  lrealed  in  [22]  and  is  much  simpler  lhan  lhe  wideband  case. 
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In  (25),  arg(-)  is  the  unwrapped  phase  operator,  which  operates  element-wise  on  a  vector.  We  note  that 
the  suboptimal  techniques  for  vibrational  MD  rely  heavily  on  peak  picking  in  the  discrete  Fourier  trans¬ 
form  of  the  returns  and  could  suffer  considerably  for  low  signal-to-noise  ratios  (SNRs). 

It  is  important  to  note  that  there  are  two  sets  of  suboptimal  estimates  for  cao,d ,  and  ^correspond¬ 
ing  to  the  two  carrier  frequencies,  whereas  only  one  suboptimal  estimate  for  Ro.  In  step  2  of  the  IRNLS, 
we  only  need  a  single  suboptimal  estimate  of  (Oa  for  initialization.  In  the  absence  of  any  a  priori  informa¬ 
tion  on  the  operating  conditions,  for  example,  the  SNR  at  the  two  carriers,  one  can  simply  average  the 
suboptimal  estimates  to  obtain  a  single  value  as  <Dosuhopt  =  (0)oUuhopl  +  Ct)o2suf)0pt)  /  2.  However,  more 

sophisticated  strategies,  such  as  Q)omhop,  =  w,0)oisuhopt  +  w2a)o2suhopn  wt+w2  =  \  can  he  used  if  prior 
knowledge  of  the  operating  conditions  is  available.  Likewise  for  d  and  <po.  It  is  further  noted  that  the 
suboptimal  estimates  for  R0  and  pi  are  not  required  to  launch  the  IRNLS,  as  these  parameters  are  not 
involved  in  step  2  maximizations  in  the  IRNLS. 


1.4.2  Rotational  MD 

Consider  the  single  blade  returns  of  (15).  Ignoring  the  noise  for  the  time  being,  the  Fourier  transform  of 
(15),  denoted  as  X^co),  is  given  by 

Xi(c»  =  //,(rw;x v')*ej4xf‘R‘lc  f)  jmJm(d)e~jm^S((0-m(0o)  (26) 


where  denotes  the  convolution  operator,  and  //.(<27,\ /)  is  the  Fourier  transform  of  /r(/?;\|/).  The  sine 
fluctuation  in  ht(n;\ j/)can  be  expressed  in  terms  of  the  spherical  Bessel  function  of  zero  order,  denoted 
by  /,(•),  leading  to  [21] 


4xfidcos(a)0n-<po) 


(27) 


\  C  J 

Expressing  jo(-)  in  terms  of  the  Bessel  function  of  the  First  kind,  and  using  its  integral  equivalent,  we 
obtain 
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h,(ir,\\i')  =  p, 


7Z 


4 nftd  cosicoji  -  (po) 


\  c 

V/r  ' 


1/2 


47rfidcos(coi>n-<po) 


(28) 


r_j  V 


r(i)ni/2)J  l. 


where  T(-)  denotes  the  gamma  function.  The  integrand  in  (28)  can  be  simplified  further  using  the  Jacobi- 


Anger  series  [21], 


cos 


i?My  i+2y(.lrJ 

c  )  °\  c 


2m 


4  nfkd 


c  J 

Using  (29)  in  (28),  and  then  applying  the  Fourier  transform,  we  obtain 


y  \cos(rn(2coon  -  2(po ))  (29) 


II  /  2 7lp  Jk  f 

Hiax,\ \>  )  = - ! - 

r(i>r(i  /  2)  J 


4  nf,d 


y  p(o)+ 


4  nftd 


y  \(S(co-2ma)it)e-2niv°  +S(co+2mcoJe2mM  ) 


2  mj<p() 


dy  (30) 


J 


U  oo 

Using  the  identity  ^Jv(x)dx  =  2^Jvt2ryl(a),  in  (30),  we  obtain 


Hi(co-,\\i')  = 


P, 


\fx  c 
X- 


r(i)r(i/2)  fd 


r= 1 

+Z(-»" 

m- 1 


S(o)- 2  mcoo)e 
+S{co+2m(t)o  )e 


~2mj<Po  \ 


2  mj<pn 


I  j. 


2ro+2r+1 ' 

r~\  C  j 


(31) 


Interestingly,  the  function  in  (31)  has  harmonics  at  0,±2coo,±4cooy.. .,  where  amplitudes  are  scaled  by  a 
complex  combination  of  Bessel  functions.  From  (26)  and  using  the  result  in  (31),  it  can  be  readily  seen 
that  the  resulting  expression  after  convolution  retains  the  original  harmonic  structure,  i.e.,  the  spectrum 
X.(cq)  has  peaks  at  0,+*y  ,±2*y  akin  to  the  vibrational  MD  spectrum.  Likewise,  it  can  be  shown 

that,  for  Q  identical  blades,  the  Fourier  spectrum  of  has  harmonics  at  0,±Qa)o,±2Qcoo,...^.  In  practice, 
similar  to  the  vibrational  MD  case,  there  exist  only  a  finite  number  of  dominant  harmonics  for  rotational 
MD.  In  essence,  the  suboptimal  estimator,  designed  for  the  vibrational  MD,  can  be  also  used  for  rotational 
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MD  to  estimate  Qa>o ,  from  which  the  suboptimal  estimate  of  the  rotational  frequency  denoted  as 
“>m.mh<v,d  =  \,2  >s  readily  obtained. 

Unlike  the  vibrational  MD  case,  it  is  not  straightforward  to  use  the  amplitude  of  the  peaks  in  the  re¬ 
turn  signal  spectrum  to  obtain  suboptimal  estimates  of  d  and  <po  for  rotational  MD.  This  is  due  to  the 
complicated  structure  of  (31)  and  (26)  involving  Bessel  functions  and  convolution.  Instead,  we  use  the 
grid  search  technique  for  estimating  d  and  <po.  The  suboptimal  estimates  for  d  and  <po  at  each  of  the 
carrier  frequencies  are  readily  obtained  as. 


> ki.u.hop, )  =  arg max  G, {&,UuhoptJ,<po) 

d  .<po 


G  t(a>„d,<pe)  = 


j  4  7tff  d  c  os(  0)on - <p0 )  Ic 


n=() 


(32) 


N- ] 


n=0 


The  cost  function  Gi(C0odJ(po)  is  easily  recognized  as  the  ML  cost  function  operating  individually  on  the 
returns  at  the  carrier  fr  Although  (32)  is  designed  for  a  single  blade,  it  can  be  used  directly  when 
multiple  blades  are  present.  In  fact,  the  presence  of  Q  blades  will  reveal  itself  as  Q  dominant  peaks  in  the 
cost  function  (32),  provided  that  Q  is  finite. 

The  grid  search  on  (32)  is  appealing  because  of  two  main  reasons.  Firstly,  as  opposed  to  traditional 
targets,  rotating  fans  do  not  assume  very  high  displacements.  Thus,  we  can  safely  assume  de  (0,1m]. 


Further,  <po  e  [0,2;r).  The  entire  parameter  space  can  then  be  evaluated  using  a  small  and  finite  grid.  For 
example,  we  can  evaluate  (32)  for  every  0.01m  in  the  d-domain,  and  similarly  for  every  0.01^  in  the  (po  - 
domain,  which  makes  the  grid  search  extremely  feasible.  Secondly,  the  visualization  of  the  cost  function 
using  the  grid  search  would  indicate  the  global  optimum  corresponding  to  one  blade,  or  several  maxima 
corresponding  to  multiple  blades  in  the  radar  returns.  In  other  words,  knowledge  of  the  number  of  blades 
can  be  readily  obtained,  if  it  is  not  known  a  priori.  For  a  reasonable  SNR,  the  cost  function  in  (32)  will 
fail  to  reveal  multiple  peaks  when  Q— >oo,  which  is  quite  unlikely  for  fans,  or  when  the  grid  size  is  too 
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large  to  accommodate  the  support  of  all  the  peaks  in  the  ( d,(po )  plane,  in  which  case  a  single  continuous 


peak  is  seen  along  the  <p„-axis. 

Note  that  the  suboptimal  techniques  for  vibrational  and  to  some  extent  rotational  MD  use  the  discrete 
Fourier  transform,  and  thus,  can  be  implemented  using  the  FFT  efficiently.  It  is  now  assumed  that  these 
suboptimal  techniques  give  rise  to  estimates  close  to  the  true  parameters  so  that  the  IRNLS-ML  iterations 
converge  to  the  optimal  ML  estimate. 

1.5  Cramer-Rao  Bounds 

The  vibrational  MD  CRBs  were  derived  in  [14].  In  this  section,  we  derive  the  multi-component  rotational 
MD  CRBs.  Let  F  be  the  Fisher  information  matrix  and  0  =  [\|/r,cr12,cr22]r  be  the  complete  parameter 
vector,  where  the  MD  parameter  vector  is  defined  as  \j /  =  [Ro,a)o,d,(po,p],p2]1  .  The  FIM  elements  can  be 

derived  in  a  compact  manner  for  the  complex  Gaussian  pdf  due  to  the  Slepian-Bangs  formula,  which  is 
given  by  [15], [25], 


(33) 


where  \i  and  C  are  the  mean  and  the  covariance  matrix  of  the  PDF,  respectively.  For  the  problem  at 
hand,  and  using  (33),  it  can  be  readily  shown  that 


(34) 


where  y .  (*|/)  and  g;<?(*j/)  have  been  defined  in  (18).  The  derivatives  with  respect  to  the  parameters  are 
required.  The  partial  derivatives  of  P(/(\y)  with  respect  to  l//r,r  =  1,2,. ..,6  are  given  by. 


(35) 
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dco„ 


•)  =  P,e 


j4xf,R0/c 


dy ,q(v')n  .  ,, ,  ,  ,.n  .  ,,n  n  .  ,  -jAnfd 

— T - o  )  +  Tj,(V  )0  gi,(V  )□  n  □  Im{^}x - 


M, 


\q  :=  \\,e‘a° , e2jm°  ]'  x e* , n  :=  [0, 1, 2. . . , N  - 1  ]' 

/ 


9  6). 


□  (V°-2) 


Re{g„(V,)}0  — — n □  Im{^  }□  ]  (36) 

c  c 

_lm(S„(v')lD  ~~~  n  □  Im(y 

,.i£MReK)1 
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where  we  define  (■)  *  as  the  element-wise  division  of  a  vector  or  a  matrix  raised  to  k- th  power,  i.e.,  if 
x  =  [x]9x2,...,xn]t ,x  “*  =[1  /  xk , l  /  xk ,...,1 1 xkNf .  Using  the  same  convention,  the  other  derivatives  which 
are  required  for  evaluating  the  FIM  are  provided  below. 
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□  (vn~2) 


M,  ,q 


P  i,(y)  =  A« 


jtxfiRJc 


fyk/iv')  n  ,  ,,  ,  ,  ,,n  ,  ,.n  .  ..  ,  j4xftd 

— 1 - 0  gi9(M')  +  Ti?(V)D  gi(? ( V  ) □  Im{^}x - i— 


<)(P„ 


\  'To 


Anftd 


Re{g,/'K  >}D  -^-Iml^in  — ^-Re{^} 

c  c 


(38) 


□  (v1-2) 


From  (33),  two  important  observations  are  in  order.  First,  the  cross  FIM  elements  with  respect  to  pt  and 
the  parameter  /?,are  zero  for  a  single  fan  blade,  and  are  in  general  non-zero  when  multiple  blades  are 


present,  i.e, 


F«*=  0,  /  =  !,  2,  <2  =  1 


(39) 


This  is  because  the  amplitude,  pt,  /  =  1,2  is  not  embedded  in  the  phase  of  the  signal,  and  similarly  the 
range  parameter  Ro  is  not  a  function  of  ht(n\\y').  The  more  important  second  observation  arises  from  the 
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fact  that  since  the  covariance  matrix  is  not  a  function  of  the  parameters  in  vg,  the  cross  FIM  elements 


with  respect  to  <j]  and  all  the  parameters  in  vy  are  zero,  i.e., 

F  v  =0,  re(l,...,6), /  =  1,2,  V(?  (40) 

Eq.  (40)  implies  that,  regardless  of  the  knowledge  (or  no  knowledge)  of  the  parameters,  cr,  the  CRBs 
depend  only  on  the  inverted  FIM  of  the  MD  parameters,  including  the  range.  Further  simplifications  of 
the  FIM  elements  are  tedious  to  derive.  The  above  equations  can  be  used  to  numerically  evaluate  the  FIM 
and  the  CRBs. 

1.6  Simulations  and  Experimental  Results 

1.6.1  Simulations 

The  carrier  frequencies  for  the  dual  frequency  operation  are  set  to  /,  =903 MHz  and  f2  =921  MHz. 
Also,  we  average  the  suboptimal  estimates  of  parameters,  such  as  coo,d,  and  (po,  corresponding  to  the 
two  carrier  frequencies,  and  use  the  averaged  values  for  IRNLS  initialization. 

We  first  consider  the  vibrational  MD.  For  simplicity  we  assume  p{  =  p.  The  SNRs  for  this  case  are 

defined  as  SNR{  =|sf||  /<J?N,\fi  =  1,2.  We  force  SNR {  =10dB  and  vary  SNR2  from  -lOdB  onwards,  in 

increments  of  5  dB.  Figures  2(a)-(d)  demonstrate  the  MSE  for  the  IRNLS  and  the  suboptimal  estimation 
schemes.  The  number  of  Monte  Carlo  trials  was  250  and  the  maximum  number  of  allowed  iterations, 
k max,  was  chosen  as  10.  The  number  of  data  samples  N  =  1024,  and  the  parameters  of  the  MD  signal 

used  were  \|/  =  [R0  =  1.3m,rf  =  0.07m, 0)o  =  0.123;r,^  =/r/3]7  The  received  signal  is  wideband  for  the 
aforementioned  carrier  frequencies.  In  the  suboptimal  estimation  scheme,  we  forced  m  =  K  =  1  which 
makes  K=[-l,0,l]r  in  eqs.  (21),  (23),  and  (25).  This  choice  of  K  provides  satisfactory  results,  and 
amounts  to  using  the  first  harmonic  and  DC  for  estimation,  i.e.,  the  suboptimal  scheme  operates  in  a 
narrowband  framework.  The  estimates  after  the  first  iteration  of  the  IRNLS  correspond  to  optimizing  the 
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2 

NLS  cost  function  ^||x-  -sl(\|i)||  /  A, and  are  also  provided  in  Fig.  2  for  comparison.  It  is  clear  from 
1=1 

Fig.  2(a)-(d)  that  the  IRNLS  gives  better  MSE  than  the  NLS  and  of  course  the  suboptimal  schemes.  The 
corresponding  CRBs  for  single  frequency  and  dual  frequency  operations  are  also  shown.  It  is  clear  that 
the  single  frequency  CRBs  guide  the  ML  behavior  for  the  dual  frequency  operations.  The  estimates  are 
clearly  above  the  dual-frequency  CRBs.  The  IRNLS  offers  superior  performance  than  the  NLS  for  all 
parameters  except  R0  for  which  both  the  NLS  and  IRNLS  give  similar  MSE.  It  is  noted  that  when 

SNRX  =SNR2  =>  <j\  =  <j\ ,  the  IRNLS  and  NLS  MSE  are  the  same  for  all  parameters.  It  is  clear  that  the 
MSE  of  the  parameters  decreases  with  the  SNR,  specifically  SNR2  as  predicted  by  the  CRBs. 

Next,  we  consider  a  rotating  fan  with  multiple  blades.  Figs.  3  and  4  shows  the  corresponding  CRBs 
with  respect  to  SNR2  and  the  data  record  length,  A,  respectively.  The  parameters  used  in  generating 

these  figures  were  \\f  =  [R0  =  1  Amyd  =  0.253my0)o  =  0.02/r,<pf>  =  K /8]7 .  The  number  of  blades  vary  from 
one  to  four.  From  Fig.  3,  it  is  evident  that  multiple  blades  yield  better  estimates  for  all  the  parameters. 
However,  from  Fig.  4,  we  observe  that,for  a  fixed  SNR  and  varying  A,  the  CRBs  corresponding  to 
multiple  blades  are  more  or  less  similar  to  those  for  a  single.  In  general,  for  the  case  of  multiple  identical 
blades,  greater  confidence  in  the  estimates  can  be  obtained  which  implies  better  estimation  of  the  parame¬ 
ters.  Consider,  for  example,  the  time-frequency  distribution  of  a  four  blade  return.  The  instantaneous 
frequency  (IF)  for  the  blades  is  identical  in  structure  but  differs  in  the  phase  parameter  cpp.  It  is  clear  that, 

for  example,  the  rotational  frequency  can  be  estimated  with  a  greater  confidence  since  four  different  IF 
curves  in  the  TF  plane  can  be  used  jointly  to  estimate  the  rotational  frequency,  rather  than  a  single  IF 
curve  when  a  single  blade  is  present,  provided  that  the  IF  curves  are  mostly  non-overlapping  in  the  TF 
plane.  The  performance  of  the  NLS  and  IRNLS  and  suboptimal  estimators  for  the  rotational  MD  corres¬ 
ponding  to  a  fan  with  three  blades  is  shown  in  Fig.  5.  Parameters  identical  to  those  for  Figs.  3  and  4  were 
used  in  Fig.  5.  For  generating  the  MSE,  the  number  of  Monte  Carlo  trials  was  set  to  200,  and  the  maxi¬ 
mum  number  of  allowed  iterations,  was  chosen  as  10.  A  coarse  grid  search  was  used  as  in  (32)  to 
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suboptimally  estimate  the  parameters  d  and  <po.  In  the  rotational  MD  simulations,  we  used  the  root- 
MUSIC  algorithm  to  estimate  the  harmonics  instead  of  peak  picking  in  order  to  compute  the  MSE  for 
SNR  values  less  than  OdB.  In  Fig.  6(a),  we  show  the  LLcost  as  given  in  (8a),  as  well  as  the  norm  of  the 
difference  in  parameter  estimates  for  various  iterations  in  the  IRNLS  in  Fig.  6(b).  Parameters  identical  to 
the  rotational  MD  simulation  of  Fig.  5  were  used  to  generate  Fig.  6,  except  that  we  show  one  Monte  Carlo 
simulation,  and  ( SNR , ,  SNR2 )  =  (1 0,  5)  dB.  It  is  clear  from  Fig.  6(b)  that  after  six  iterations  of  the  IRNLS, 

the  norm  becomes  zero.  In  essence,  there  exists  a  cluster  or  accumulation  point  as  shown  in  Fig.  (6b).  The 
discussion  in  Appendix-A  gives  more  details  on  the  cluster  point  with  respect  to  the  IRNLS. 

1.6.2  Experimental  Results 

In  this  section,  we  present  results  of  IRNLS  applied  to  various  MD  signal  returns  measured  using  a  dual¬ 
frequency  radar  in  a  laboratory  setting.  First,  a  12in  conducting  sphere  was  tied  to  the  ceiling  with  a  nylon 
string,  and  excited  to  oscillate  in  a  simple  harmonic  fashion.  This  experiment  corresponds  to  the  vibra¬ 
tional  MD  scenario.  The  approximate  range  to  the  target  was  16.4ft  or  4.99m,  the  carrier  frequencies  were 
chosen  to  be  906.3  and  919.6  MHz,  and  the  sampling  frequency  was  100Hz.  More  details  on  the  experi¬ 
ment  can  be  found  in  [26].  We  processed  only  the  First  seven  seconds  of  data,  comprising  700  samples,  in 
order  to  avoid  damped  oscillations.  Figure  7  shows  the  IRNLS,  NLS,  and  suboptimal  instantaneous 
frequency  (IF)  trajectories  overlaid  on  the  spectrogram  of  the  data.  Clearly,  the  IRNLS  yields  better 
estimates  when  compared  to  the  NLS,  and  agrees  with  the  IF  of  the  data  for  both  carrier  frequencies.  The 
range  estimates  for  IRNLS  and  NLS  were  4.89m  and  4.83m,  respectively. 

Second,  the  IRNLS  is  applied  to  measured  returns  from  a  table-top  fan  with  4-metallic  blades.  The 
carrier  frequencies  are  chosen  as  903.6  MHz  and  921  MHz.  The  distance  from  the  center  of  the  fan  to  the 
radar  was  3  m,  and  the  sampling  frequency  was  chosen  to  be  5  kHz  to  avoid  aliasing.  Two  datasets,  one 
for  azimuth  equal  to  0°,  and  the  other  for  60°  azimuth,  were  used  for  the  ML  analysis.  Fig.  8(a)  shows  the 
spectrogram  of  the  raw  data  for  0°  azimuthal  aspect,  using  a  rectangular  window  of  101  samples  which 
corresponds  to  around  130  %  of  the  period  of  the  rotation.  For  subsequent  ML  analysis,  the  data  was 
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decimated  by  a  factor  of  5.  The  2D  grid  search  cost  function,  as  given  in  (32)  was  used,  and  the  cost 
function  is  shown  in  Fig.  8(b)  for  the  second  carrier.  The  grid  was  initialized  with  the  rotational  frequency 
estimate  obtained  from  the  Fourier  spectrum.  We  can  clearly  see  4  dominant  peaks  corresponding  to  the  4 
blades.  The  suboptimal  estimates  were  used  to  initialize  the  ML  algorithm.  The  range  estimates  and  final 
LL  costs  are  provided  in  Table  I.  The  results  in  Table  I  indicate  that  both  the  NLS  and  IRNLS  perform 
exceptionally  well  for  this  dataset.  Fig.  9(a)  shows  the  spectrogram  of  the  raw  data  (not  resampled), 
corresponding  to  60°  azimuthal  aspect,  using  the  same  window  length  as  in  Fig.  8(a).  The  2D  grid  search 
cost  function  corresponding  to  the  second  carrier  frequency  is  shown  in  Fig.  9(b).  Similar  to  Fig.  8(b),  the 
2D  grid  in  Fig.  9(b)  clearly  shows  4  dominant  peaks.  The  corresponding  ML  estimates  for  the  range  are 
shown  in  Table  II,  which  shows  that  the  IRNLS  performs  better  than  the  NLS. 

1.7  Conclusions 

In  this  chapter,  we  considered  dual  frequency  Doppler  radars  for  range  estimation  of  moving  targets  with 
application  to  urban  sensing.  The  ML  estimator  was  derived  for  the  micro-Doppler  motion  profile,  which 
is  commonly  exhibited  by  indoor  moving  targets.  It  was  shown  that  the  ML  estimator,  although  not 
solvable  in  closed  form,  can  be  provided  using  a  step-wise  iterative  algorithm  termed  as  the  IRNLS.  The 
algorithms  solution  and  procedure  depends  on  the  motion  profile  model.  The  initial  conditions  are 
provided  using  suboptimal  values  which  rely  on  the  harmonic  nature  of  the  radar  returns.  For  simulated 
data,  the  algorithm  was  shown  to  be  superior  in  terms  of  the  MSE  when  compared  to  suboptimal  estima¬ 
tors.  The  iterative  ML  was  also  applied  to  data  measurements  corresponding  to  indoor  moving  targets, 
yielding  superior  IF  estimates  and  good  estimates  of  the  range. 
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Appendix-A:  Convergence  of  the  IRNLS  algorithm 

In  this  appendix,  we  do  not  assume  any  specific  motion  profile.  The  results  are  general  and  applicable  to 


0  =  [\gr,(j,2,<J2  ]7,vj;  =  [^p^2’^3 . VPV  for  the  IRNLS. 


Proposition  1:  The  FIM  for  the  dual  frequency  radar  can  be  expressed  as 


F(V)  0  0  1 


F(0)  =  0  N  /  <j|4  0 

0  0  N I  o  I  \ 


(A.l ) 


Proof:  The  proof  of  the  above  proposition  is  standard  and  can  be  easily  obtained  using  the  Slepian-Bangs 


formula  [1 5], [25]. 

Proposition  2:  If  the  minimization  in  (10)  is  carried  out  using  the  method  of  Fisher  scoring,  with  step  size 
equal  to  unity,  then  the  IRNLS  algorithm  becomes  the  Fisher  scoring  technique  for  the  ML  problem. 
Proof:  We  first  note  that  for  a  fixed  arbitrary  covariance  matrix,  the  gradient  vector  with  respect  to  the 
desired  parameter  vector  \p  of  the  negative  LL  is  identical  to  the  gradient  vector  of  the  cost  function  in 

(10).  Hence,  the  Fisher  scoring  technique  for  the  ktf  iteration  for  estimating  \j t  k  is  then  given  as  [27] 


(A. 2) 


where  F  '(V*  i)  and  )  are  the  FIM,  and  the  gradient  vector  of  the  cost  function  in  (10), 

obtained  using  the  previous  estimate,  and  Xk  is  the  step  size  parameter.  The  negative  sign  in  (A. 2) 
preceding  the  second  term  indicates  the  minimization  of  the  cost  function  (10),  or  equivalently  minimiza¬ 
tion  of  the  negative  LL  .  The  Fisher  scoring  applied  to  the  negative  LL  for  the  entire  parameter  vector  0, 
takes  the  following  form. 


o*=o*-,-A 


o 

0 


0 


o  lr  -vj(V>  i 

0  —3  ln(/?(x;\|/)  /  Sex,2 


(A. 3) 


0 


0-20-1) 7  N  jL-9  ln(/?(x;  V)  /  9cr2 1 


2 


The  sign  is  Hipped  for  the  conventional  Fisher  scoring  technique  which  is  applied  to  the  maximizing  the  LL  [27]. 
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The  term  -9ln(p(x;i|/))/ dcrf ,  V/  =  1,2  was  an  intermediate  step  in  deriving  (7).  Using  this  intermediate 
step  with  Ak  =  1,  and  multiplying  the  matrix  and  the  gradient  vector  in  (A3),  we  obtain  (A. 2)  with 
Xk  =1,  and 

=|lx<  -s.('t'*-i)||2/jV’  V'  =  *>2  (A.4) 

Equation  (A.4)  is  the  remaining  key  equation  in  the  IRNLS  for  updating  in  (10).  Q.  E.  D. 

The  special  case  of  Xk  =1  is  oftentimes  used  in  the  literature  (see  [15],  and  also  [27]).The  above 
proposition  implies  that  the  IRNLS  is  a  special  case  of  Fisher  scoring  algorithm.  However,  as  well 
known,  the  Fisher  scoring  converges  if  the  suboptimal  estimates  used  to  initialize  it  are  close  to  the 
optimal  ML  solution.  In  [27,  remark  4.1  and  remark  4.2],  under  certain  conditions  it  was  shown  that  the 
Fisher  scoring  method  converges  to  at  least  a  local  minimum  of  the  negative  LL  as  N  — >  ,  i.e.,  asymp¬ 
totically.  The  same  follows  from  proposition  2  for  the  IRNLS,  more  details  can  be  seen  in  [27].  If  the 
negative  LL  objective  is  convex,  which  is  a  standard  assumption  for  ML  techniques,  then  the  Fisher 
scoring  converges  surely  to  the  global  minimum  [28].  Convergence  is  also  attained  if  the  initialization  is 
contained  within  a  set  including  the  optimal  solution,  where  the  objective  is  locally  convex  [28].  For  step¬ 
wise  iterative  ML  schemes,  an  interesting  result  is  provided  in  [29,  lemma-1],  where  the  authors  prove  the 
existence  of  a  cluster  point  or  an  accumulation  point,  denoted  as  O',  for  the  sequence  of  parameter 
iterates  given  by  {0A  }  if  the  negative  LL  is  bounded  below,  which  is  definitely  satisfied  by  the  Gaussian 
pdf.  We  can,  thus,  deduce  that  the  IRNLS  iterates  also  exhibit  similar  behavior.  Essentially,  if  there  exists 
a  unique  cluster  point  O'  for  any  sequence  {0*},then  0'  is  the  ML  solution  0(x).  This  directly  implies 

that  hm^^lO^,  -0*1  =  0.  The  existence  of  the  cluster  point  was  also  proven  in  [28]  for  the  constrained 

Fisher  scoring  technique.  We  now  prove  that  the  Fisher  scoring  and  hence  the  IRNLS  increases  the  LL  for 
every  iteration  monotonically. 

Proposition  3:  Under  certain  assumptions,  the  Fisher  scoring  and  hence  the  IRNLS  increases  the  LL  for 

every  k ,  until  convergence  as  N  — 
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Proof:  It  is  required  to  prove  1n(p(x;0t+l))-ln(p(x;0*))>O.  Expanding  the  pdf  about  the  arbitrary  point 


AkF  '(0t)VJ(0l)t  where  VJfO^)  is  the  gradient  of  the  LL  evaluated  at  the  previous  0*,  using  the 
Taylor  series  expansion. 


ln(x;0t+1  )-ln(x;0t )  =  Ak  VrJ(6*  )F  <0*  )VJ(0* ) 

A2 

+  yVrJfOt  )(Fr(0*  ))"1[V(VrJ(6Jt  ))]VJ(0*  )F*'  (0, ) 


(A. 5) 


where  [V(V' J(0t ))]  is  the  Hessian  matrix  evaluated  at  0*.  Now  using  the  results  in  [30],  i.e.,  as 
N  — >  °o  the  negative  Hessian  matrix  converges  to  the  FIM  with  probability  1,  we  can  replace  the  negative 
Hessian  matrix  with  the  FIM,  F(0*)  in  (A. 5),  further  assuming  that  the  LL  can  be  approximated  as  a 
quadratic,  i.e.,  ignoring  the  higher  order  terms,  we  get 


A2 

ln(x;0*+1  )-ln(x;0* )  =  Ak Vr  J(0*  )F  1  (0t  )VJ(0*  Vr  J(0t  )F  1  (0*  )VJ(0; ) 

=  (Ak  -A2k  / 2)Vr  J(6t  )F-' (0k )VJ(0t ) 


Since  the  FIM  is  nonsingular  and  positive  definite  V0*,  with  /^e(0,l],  we  have  (A. 6)  is  always 
positive,  and  becomes  zero  at  the  local  optimum  solution.  This  proves  the  proposition  for  the  Fisher 
scoring  approach.  Likewise,  substituting  Ak  =1  ,\/k  it  is  readily  seen  that  (A. 6)  is  always  positive.  This 
proves  the  proposition  for  the  IRNLS.  Q.  E.D. 

In  this  Appendix,  we  have  proved  that  the  IRNLS  is  a  special  case  of  Fisher  scoring.  Therefore, 
some  of  the  results  for  Fisher  scoring  directly  follow  through  to  the  IRNLS.  In  the  simulations,  we 
demonstrate  the  existence  of  the  cluster  point  and  the  monotonic  increase  in  the  LL  objective  for  every 
iteration. 
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Fig.  1  Flowchart  of  IRNLS  ML. 


(d) 
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(c) 


Fig.  2  MSE  of  vibrational  MD  parameters.  IRNLS  Comparsion  wrt  CRBs  and  NLS. 
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Fig.  3  CRBs  for  multiple  blades,  N=5\2,  SNR|=15dB,  x-axis  SNRi  dBs,  y-axis  CRBs,  blue  (  1  blade),  red 

(2  blades),  black  (3  blades),  magenta  (4  blades). 
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Fig.  4  CRBs  for  multiple  blades,  SNR|=15dB,  SNR2  =-I0dB,  x-a\is  N,  data  record  length,  y-axis  CRBs, 
blue  (  1  blade),  red  (2  blades),  black  (3  blades),  magenta  (4  blades). 


Fig.  5  MSE  performance  of  rotational  MD  parameters,  SNRi=5dB,  SNR2  is  varied. 
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Fig.  6  Figures  showing  the  (a)  LL  cost,  (b)  cluster  point. 


.  Time  m  $ 


1  2  3  4  5  6 


@1^.  Time  in  s 

Fig.  7  Spectrogram  of  measured  radar  returns  for  vibrational  MD,  along  with  the  final  IF  fits. 
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Fig.  8(a)  Spectrograms  of  the  measured  returns  at  0°  azimuth  at  the  two  carriers. 
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Fig.  8(b)  LL  grid  (0°  azimuth  )  using  single  blade  cost  function,  for  carrier  f2.  Four  blades  are  clearly 

seen.  Number  of  samples  used  is  1200. 

Table  I.  Range  estimates  at  0°  azimuth.  True  range  is  3m. 


Data  chunks 

Minimized 

Minimized 

Range  estimate 

Range  esti- 
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comprising  of 

LLcost  (NLS) 

LLcost  (IRNLS) 

(NLS)  in  meters 

mate  (IRNLS) 

500  samples. 

from  8(a) 

from  8(a) 

in  meters 

100th  sample 

-28.7104 

-28.9614 

3.0046 

2.9977 

onwards 

700'"  Sample 

-28.5853 

-28.9000 

3.0052 

2.9985 

onwards 

1400th  Sample 

-29.6069 

-30.0953 

3.0132 

3.0051 

onwards 

Table  II.  Range  estimates  at  60°  azimuth.  True  range  is  3m. 


Data  chunks 

comprising  of 

500  samples. 

Minimized 

LLcost  (NLS) 

from  (8a) 

Minimized 

LLcost  (IRNLS) 

from  8(a) 

Range  estimate 

(NLS)  in  meters 

Range  esti¬ 
mate  (IRNLS) 

in  meters 

100th  sample 

onwards 

-31.8872 

-31.9728 

3.2449 

3.2336 

700th  Sample 

onwards 

-31.9589 

-32.1712 

3.1552 

2.9286 

1400th  Sample 

onwards 

-31.7151 

-31.7550 

3.0900 

3.087 
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Fig.  9(a)  Spectrograms  of  the  measured  returns  at  60°  azimuth  at  the  two  carriers. 


Fig.  9(b)  LL  grid  (60°  azimuth)  using  single  blade  cost  function,  for  carrier  f2.  Four  blades  are  clearly 

seen.  Number  of  samples  used  is  1200. 
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2.  Dual  Frequency  Doppler  Radars  for  Indoor  Range  Estimation:  Cramer- 

Rao  Bound  Analysis 


Abstract 

Single  frequency  Doppler  radars  are  effective  in  distinguishing  moving  targets  from  stationary  targets,  but 
suffer  from  inherent  range  ambiguity.  With  a  dual-frequency  operation,  a  second  carrier  frequency  is 
utilized  to  overcome  the  range  ambiguity  problem.  In  urban  sensing  applications,  the  dual-frequency 
approach  offers  the  benefit  of  reduced  complexity,  fast  computation  time,  and  real  time  target  tracking. 
We  consider  a  single  moving  target  with  three  commonly  exhibited  indoor  motion  profiles,  namely, 
constant  velocity  motion,  accelerating  target  motion,  and  micro-Doppler  motion.  RF  signatures  of  indoor 
inanimate  objects,  such  as  fans,  vibrating  machineries,  and  clock  pendulums,  are  characterized  by  micro- 
Doppler  motion,  whereas  animate  translation  movements  produce  linear  FM  Doppler.  In  this  chapter,  we 
derive  Cramer-Rao  bounds  (CRB)  for  the  parameters  defining  indoor  target  motions  under  dual-frequency 
implementations.  Experimental  data  is  used  to  estimate  micro-Doppler  parameters  and  to  validate  the 
CRBs. 
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2.1  Introduction 


The  emerging  area  of  through-the-wall  sensing  addresses  the  desire  to  see  inside  enclosed  structures  and 
behind  walls  in  order  to  determine  building  layouts,  scene  contents,  and  occupant  locations.  This  capabili¬ 
ty  has  merits  in  a  variety  of  civilian  and  military  applications,  as  it  provides  vision  into  otherwise 
obscured  areas,  thereby,  facilitating  information-gathering  and  intelligent  decision-making  [1-4]. 

Motion  detection  is  a  highly  desirable  feature  in  many  through-the-wall  applications,  such  as  fire¬ 
fighters  searching  for  survivors  or  law-enforcement  officers  involved  in  hostage  rescue  missions.  Doppler 
discrimination  of  movement  from  stationary  background  clutter  can  be  readily  used  in  these  applications. 
Such  systems  can  be  classified  into  zero,  one,  two,  or  three  dimensional  systems  [1].  A  0-D  system 
simply  provides  a  go/no-go  motion  detection  capability  and,  as  such,  can  employ  single  frequency 
continuous  wave  (CW)  waveforms  for  scene  interrogation.  A  1-D  system  can  provide  range  to  moving 
targets  by  employing  modulated  or  pulsed  signals.  A  2-D  system  provides  both  target  range  and  azimuth 
information,  whereas  a  3-D  system  adds  the  dimension  of  height  to  the  offerings  of  a  2-D  system. 
However,  the  higher  level  of  situational  awareness  offered  by  2-D  and  3-D  systems  is  obtained  at  the 
expense  of  increased  hardware  and  software  complexity  and  higher  computational  load.  A  1-D  system 
provides  a  good  compromise  between  level  of  situational  awareness  and  cost  versus  size  and  weight.  This 
is  specifically  the  case  when  localizing  moving  targets  is  of  primary  interest. 

In  this  chapter,  we  consider  a  1-D  system  for  range-to-motion  estimation  based  on  Doppler  radars. 
Single  frequency  CW  radars  suffer  from  range  ambiguity  [5],  This  problem  is  more  pronounced  for 
through-the-wall  applications  since  the  carrier  frequencies  are  in  the  few  GHz  range  due  to  antenna  size 
and  frequency  allocation  management  issues.  An  additional  carrier  frequency  can  be  introduced  to  resolve 
the  uncertainty  in  range.  The  two  carrier  frequencies  can  be  chosen  such  that  the  maximum  unambiguous 
range  is  either  equal  to  or  greater  than  the  spatial  extent  of  the  urban  structure.  In  this  case,  the  true 
solution  is  the  one  which  corresponds  to  the  target  inside  the  enclosed  structure.  It  is  noted  that  other 
radar  techniques,  such  as  the  swept  frequency,  pulse  compression  etc.,  can  be  used  to  reduce  the  uncer¬ 
tainty  in  range  [5].  However,  the  operational  logistics  and  system  requirements  for  through-the-wall 
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operations,  such  as  cost,  hardware  complexity,  and  portability,  may  prohibit  the  use  of  such  techniques. 
The  dual-frequency  approach,  although  not  new,  meets  all  of  the  above  requirements  and  is  likely  to 
emerge  as  a  leading  forerunner  in  modem  urban  sensing  systems  [6-9]. 

We  focus  on  three  target  motion  and  range  profiles,  namely,  linear  translation,  acceleration,  and  the 
micro-Doppler  (MD)  motions  that  are  commonly  encountered  in  urban  sensing  applications.  For  example, 
the  swinging  of  the  arms  and  the  legs  of  a  person  during  walking  can  be  modeled  as  micro-Doppler, 
whereas  the  gross  motion  of  the  torso  can  be  represented  by  translation.  We  develop  Cramdr-Rao  bounds 
(CRBs)  for  the  parameters  corresponding  to  the  aforementioned  three  indoor  motion  profiles.  We  provide 
the  CRB  for  the  range  estimate  in  addition  to  motion  speed,  harmonic  frequency,  and  chirp  parameters. 
Specifically,  we  show  that  the  dual  frequency  approach  increases  the  Fisher  information  compared  to  a 
single  frequency  operation  and,  therefore,  lowers  the  CRBs  [10].  In  the  analysis,  the  radar  aspect  angle  6 
is  incorporated  into  the  general  CRB  derivations.  When  the  aspect  angle  is  assumed  unknown,  the  FIM 
becomes  singular,  rendering  the  parameters  non-identifiable.  It  is  noted  that  the  main  contribution  of  the 
chapter  is  the  derivation  of  the  CRBs,  and  no  estimation  technique  is  advocated  here.  For  MD  motion,  we 
use  experimental  data  to  estimate  micro-Doppler  parameters  and  to  validate  the  CRBs.  We  employ 
suboptimal  techniques  that  make  use  of  the  inherent  impulsive  structure  of  the  spectrum  of  the  returns  to 
estimate  some  of  the  parameters  in  the  MD  model  [10].  Maximum  Likelihood  optimal  estimator  for  MD 
motion  profile  is  presented  in  [10], 

A  brief  outline  of  the  chapter  is  as  follows.  In  Section  2.2,  the  signal  model  is  introduced  along  with 
the  governing  range  equations  for  the  three  motion  profiles.  In  Section  2.3,  a  generic  derivation  of  the 
CRB  is  presented,  incorporating  the  nuisance  parameters,  and  subsequently  the  CRBs  of  the  three  motion 
profiles  are  derived.  Supporting  simulation  and  experimental  results  are  provided  in  Section  2.4.  Section 
2.5  contains  the  conclusions. 
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2.2  General  Signal  Model 

Below,  we  briefly  discuss  the  range  ambiguity  problem,  and  introduce  the  motivation  for  adopting  carrier 
frequency  diversity  in  CW  radar  for  range  estimation. 


2.2.1  Range  ambiguity  in  dual-frequency  radar 

For  a  dual  frequency  radar  employing  two  known  carrier  frequencies  fx  and  /2,  the  baseband  signal 
returns  due  to  a  single  moving  target  are  expressed  as, 

(P\(nTs)  =  4nf\R(nTs)l  c  ,  <p2(nTs)  =  4x f2R(nTs) /  c  (]) 

where  Ts  is  the  sampling  period,  R(nTs)  is  the  range  of  the  target  at  sample  n ,  and  c  is  the  speed  of 
light  in  free  space.  Dropping  Ts  in  the  phase  and  range  notations  for  convenience,  and  without  consider¬ 
ing  phase  wrapping,  the  range  of  the  target  is  given  by  [6] 

KQ.cIftW-ftW)  (2) 

4*(/2-/,) 

In  reality,  however,  phase  observations  are  wrapped  within  the  [0,2;r)  range.  Therefore,  the  true  phase 
can  be  expressed  as 


0{!rue)  (n)  =  (p2  ( n )  -  ( n )  +  Imir ,  (3a) 

where  m  is  an  unknown  integer.  Accordingly,  the  range  estimate  is  subject  to  range  ambiguity  [7],  i.e, 

(3b) 


R(n)  =  c[</>2  (n)  -  (az)]  +  cm 


W2-/,)  2(/2  -  /j) 

The  second  term  in  the  above  equation  induces  ambiguity  in  range.  For  the  same  phase  difference,  the 
range  can  assume  infinite  values  separated  by 

/?w=c/(2(/2-/,))  (3c) 

which  is  referred  to  as  the  maximum  unambiguous  range.  From  the  above  equation,  Ru  is  inversely 
proportional  to  the  difference  of  the  carrier  frequencies.  Consider  a  single  frequency  Doppler  radar 
operating  at  1  GHz.  In  this  case,  Ru  -  c!2fx  =15  cm.  For  dual  frequency  operation  at  /,  =IGHz, 
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and  f2  =  1.01  GHz,  Ru  =  15  m,  according  to  (3a).  The  increase  in  Ru  from  15  cm  to  15  m  is  sufficient  for 
target  location  in  rooms  and  small  buildings.  The  other  multiple  range  solutions  given  by  (2),  and  (3b)  are 
simply  rejected  since  they  do  not  fall  within  the  spatial  extent  of  the  urban  structure.  In  the  following 
sections,  we,  therefore,  assume  m  =  0,  for  simplicity.  By  choosing  closely  separated  carrier  frequencies, 
any  maximum  unambiguous  range  can  be  achieved.  Coherent  phase  estimation  may,  however,  be  com¬ 
promised  if  the  frequency  difference  is  too  small  to  overcome  noise  effects  and  frequency  drifts  in  down 
conversions  [11].  Performance  analysis  of  the  dual-frequency  radar  in  the  presence  of  noise  and  frequen¬ 
cy  drifts  is  provided  in  [  1 2]-[  1 3]. 

It  is  noted  that  in  the  presence  of  multiple  moving  targets,  the  phase  extracted  from  the  radar  returns 
cannot  be  used  for  range  estimation  of  each  target  separately.  This  is  because  the  phase  terms  correspond¬ 
ing  to  different  targets  are  superimposed  and  cannot  be  separated  in  the  time  domain.  This  drawback  of 
the  dual-frequency  scheme  can  be  overcome  provided  that  the  target  Doppler  signatures  are  separable  in 
the  frequency  or  the  time-frequency  (TF)  domain  [14].  With  the  extraction  of  the  phase  of  the  individual 
targets  in  the  Doppler  domain,  (2)  can  then  be  used  to  obtain  the  corresponding  target  range.  Likewise  for 
multipath  which  is  similar  to  multiple  targets  the  same  procedure  can  be  applied,  provided  they  are 
separable  in  the  frequency  or  time-frequency  domain. 


2.2.2  Signal  model  for  CRB  derivation 

For  a  dual  frequency  radar  employing  two  known  carrier  frequencies  /,  and  /2,  the  baseband  signal 
returns  due  to  a  single  moving  target  are  expressed  as. 


f 

x,  ( n )  =  xx  (nTs  )  =  (/?,+  /i,  (/?; y ))  exp 

\ 

x2(n)  =  x2(nTs)  =  (p2+h2(n\\Y))zxp 
{«|/ie  Z+,/j  =  0,1, 2,3. ../V - 1} 


An  Win)' 
c  j 
'  4nf2R(n) 

s  c 


+  V,(/2)  =  J,(/2)  +  V1(/2) 

\ 

+  v2(n)  =  s2(n)  +  v2(n) 

J 


(4) 


where  pt  is  the  amplitude  of  the  return  at  frequency  fi,i  =  \,2.  The  scalar  function  dependent 

on  both  the  sample  index  n  and  a  vector  y  of  desired  parameters,  is  added  to  represent  possible  cyclic 
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amplitude  fluctuations  for  the  MD  model.  It  is  noted  that  for  the  case  of  constant  Doppler  and  accelerating 
target  motion  profiles,  ht{n\ vj/)  =  0,  since  for  these  targets  the  radar  cross  section  (RCS)  can  he  assumed 

constant  for  a  relatively  small  time  on  target.  Z+  is  the  domain  of  positive  integers.  The  two  noise 
sequences,  corresponding  to  the  two  frequencies  of  operation,  over  the  observation  period  N,  are 

complex  additive  white  Gaussian  noise  (AWGN)  and  uncorrelated,  i.e.,  vi(n)  =  v](nTs)  -  CN(0frr^ ), 

v2(n)  =  v2{nTs)  ~  CN( 0,0%  ),  E{v{(n)v2(n)}  =0,  where  <T,Z  and  o\  represent  the  noise  variance  for  the 
two  frequencies,  respectively  and  is  the  complex  conjugate  operator.  To  aid  in  the  derivation,  we 
vectorize  the  observations  and  describe  the  signal  returns  by  a  joint  probability  density  function  (pdf).  The 
received  signals  at  the  two  frequencies  are  appended  to  form  a  long  vector, 

x  =  s  +  v  =  [x1(0),*,(l),*l(2)....x1(N-l),*2(0),.x2(l) . *2(yV-l)f  =[x,  x2]'  (5) 

where 

S  =  [5,  (0),  5,  (1),  5,  (2)... .5,  (N  -  1),  S2  (0),  52 (1 ) S2(/V-l)]r  =[s,  s2 f 

(6) 

v  =  [ V,  (0),  V,  (1),  V, (2)....v,  (N- 1),  v2 (0),  v2 (l)...v2 (N  -  l)f  =  [ v,  v2f 
The  mean  and  covariance  matrix  of  x  are  given  by 


£{x}  =  H  =  [s,  s2]\  £{(x-p)(x-p)w }  =  C  = 


JNxN 


all 


(7) 


The  joint  pdf  of  the  data,  assuming  a  multivariate  complex  Gaussian  distribution,  is 


£(x;s)  =  2 ,v  T  exp(-(x  -  ft)"  C  1  (x  -  n)>  (8) 

K  '  det(C) 

1)  Constant  Doppler  Frequency 

For  the  constant  Doppler  frequency  model,  the  target  range  is  parameterized  as, 

R(n)  =  Ro  +  vcos(d)n,  ^  (9) 


where  Ro  is  the  initial  range  of  the  target  at  time  n  =0.  The  discrete  velocity  v  of  the  target  is  its  actual 
velocity  multiplied  by  the  sampling  period.  The  parameter  6  is  the  viewing,  or  aspect,  angle  of  the  radar 
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with  respect  to  the  target.  In  other  words,  the  component  vcos(#)  is  the  identifiable  Doppler  component 
registered  by  the  radar  system.  Inserting  (9)  into  (4),  it  is  clear  that  the  constant  Doppler  model  is  analog¬ 
ous  to  the  sinusoidal  parameter  estimation  [15,  16].  Hence,  the  identifiability  criteria  of  the  parameters  are 
similar  to  the  traditional  single  frequency  estimation  problem.  However,  some  details  must  be  added.  The 
identifiability  criteria  are  developed  by  putting  0~Q  in  (9)  because  for  any  other  value  of  6  in  the 
defined  range,  the  criteria  are  valid.  The  terms  2/,v/c  and  2f2v/c  in  the  phase  of  the  signal  returns  are 
the  Doppler  frequencies,  analogous  to  the  frequencies  (radian/s)  in  the  sinusoidal  parameter  estimation 
problem.  Hence,  for  v  to  be  identifiable,  -0.5  <  2/jv/e  <  0.5  and  -0.5  <  2/2v/c  <  0.5  [15].  Since  we 
have  assumed  /2  > /}  implicitly,  -c/4/2  <v<c/4f2  represents  the  bounds  on  the  positive  and  negative 
Doppler,  respectively.  For  example,  if  f2  =  1  GHz,  then  the  discrete  velocity  is  constrained  to  a  maximum 

of  0.075,  assuming  a  nominal  sampling  frequency  of  100Hz  and  the  analog  velocity  is  7.5  m/s.  It  is  noted 
that  indoor  animate  moving  targets  such  as  humans  and  pets,  do  not  assume  such  high  velocities.  Typical¬ 
ly,  in  indoor  settings,  we  expect  analog  velocities  below  lm/s  which  allows  the  sampling  frequency  to  be 
reduced  to  approximately  15  Hz.  The  terms  4 7rf]R()/c  and  4 7rf2R(Jc  represent  the  phase  in  the  sinu¬ 
soidal  parameter  estimation  model.  Clearly,  Ra  can  be  uniquely  estimated  if  the  maximum  unambiguous 
range  is  larger  than  the  building  extent,  as  determined  by  (3c).  For  this  motion  profile,  as  mentioned 
earlier,  we  use  /zf  (zz; \(/)  =  0,  since  the  RCS  is  assumed  to  be  constant. 

2)  Micro-Doppler  (MD) 

The  range  for  this  motion  profile  is  parameterized  as  [17], 

R(n)  =  Ro+dcos(0)cos{co<tn-<po),  (d,coo)*  0,  e  e  0°) 

In  the  above  equation,  the  parameter  COo  is  the  discrete  version  of  the  rotational  or  the  vibrational  fre¬ 
quency.  The  parameter  dcos(6)  is  the  maximum  displacement  from  the  initial  range  R0  along  the  radar’s 
line  of  sight  (LOS).  For  the  specific  case  of  rotation,  2dcos(6?)  represents  the  diameter  of  the  circular 
trajectory  in  the  direction  of  the  LOS.  The  parameter  <po  denotes  the  initial  phase  of  the  MD  motion.  The 
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MD  signal  returns  represent  the  traditional  sinusoidal  FM  model  [10,  11],  and  are  obtained  by  substituting 
(10)  in  (4).  It  is  noted  that  the  sinusoidal  FM  model  is  a  good  fit  for  motion  profiles  of  indoor  inanimate 
objects,  such  as  fans,  vibrating  machineries,  and  clock  pendulums.  For  the  more  complicated  human 
motion,  the  sinusoidal  FM  model  is  a  rough  approximation  to  the  swinging  of  the  arms  and  the  legs 
during  walking. 

The  MD  model  can  be  further  decomposed  into  vibration  and  rotation  models  based  on  the  scalar 
function  ht(n\\\i).  For  typical  indoor  targets  undergoing  pure  vibration  (i.e.  the  target  moves  back  and 
forth),  the  displacements  are  small  relative  to  the  target  range,  especially  for  longer  radar  standoff 
distances  from  the  wall,  and  the  RCS  is  considered  to  be  constant,  i.e.,  ht(n;\\i)  =  0.  On  the  other  hand, 
for  a  rotating  target,  a  fixed-location  radar  observes  different  elevation  aspects  of  the  target,  thereby 
inducing  a  cyclic  amplitude  modulation  in  the  radar  returns.  These  amplitude  fluctuations  are  indeed 
dependent  on  the  target  geometry.  In  this  case,  the  scalar  function  h^rr^)  can  be  modeled  as 

hl(n\\\>)  =  A Pi  cos  (coon  -<pt)  (11) 

where  the  parameter  A p{  represents  the  maximum  deviation  from  the  nominal  amplitude  piy  whereas  the 
parameter  (j)t  e  [0,  2tt)  represents  the  phase  at  time  n  =0.  As  an  example,  consider  a  rotating  fan  with  one 
blade,  at  one  extreme  instant  of  time  the  tip  of  the  blade  is  seen  representing  a  near  zero  RCS,  whereas  at 
another  extreme  time  instant,  the  entire  length  of  the  blade  is  seen  by  the  radar  representing  a  maximum 
RCS.  The  cyclic  amplitude  modulations  will  be  demonstrated  experimentally  in  the  simulations  section 
and  are  in  agreement  with  the  model  in  (1  1 ). 

The  identifiability  criteria  for  MD  signals  are  developed  by  substituting  0  =  0  in  (10).  For  the  varia¬ 
ble  RoJ  this  case  is  identical  to  the  case  of  constant  Doppler.  However,  it  is  not  straightforward  to  derive 
the  parameter  identifiability  criteria  for  the  rest  of  the  parameters  using  the  time  domain  model  alone,  as 
was  the  case  in  the  constant  Doppler  model.  We  use  the  Fourier  transform  to  aid  in  the  derivation  and  to 
shed  more  light  on  the  identifiability  bounds  and  the  discrete  parameterization  in  the  MD  model.  The 
noise  free  MD  returns  are  given  by, 
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si(n)  =  (pi+hs(n;\ |/))exp 


•  4/efjR0  .  4nftd cos(coon  -  <po ) 


■+y 


Vi  =  1,2 


(12) 


Let  Jk(-)  be  the  klh  order  Bessel  function  of  the  first  kind.  Accordingly  [18,  19], 


exp|  j^f-cos(a>0n-<po)  ]=  /  exp  (jk(a>0n-<p0))Jk 


*=- 


4  nf.tA 

\  c  ) 


(13) 


Using  (12)  and  the  Fourier  representation  of  (13),  and  ignoring  ht(n\\\f),  as  it  does  not  affect  identifiabil- 
ity,  we  obtain 


S'(eJ  )  =  p,  exp 


j4*fiK 

c 


X  )kS{(o-kco0)t\^>{-jk(Po)Jk 


k=- 


4k  \ 


,  Vi  =  1,2 


(14) 


h 

V  c  ) 

As  in  traditional  multiple  sinusoidal  parameter  estimation,  for  the  frequencies  to  be  identifiable,  they  must 
lie  in  the  interval  [-/r,/r).  This  condition,  which  translates  in  the  underlying  case  to 

~K<k(Oa  < K ,  V/:,  ke  Z,  cannot  be  imposed,  since  there  is  an  infinite  number  of  harmonics.  Moreover, 
it  is  clear  from  (14)  that  the  signal  is  not  analytic  due  to  the  presence  of  harmonics  in  the  negative 
frequencies.  However,  it  is  well  known  that  the  Bessel  functions  decrease  rapidly  with  k.  Therefore,  we 
impose  the  following  two  assumptions  to  derive  the  identifiability  constraints:  1)  Most  of  the  energy  in 
the  signal  can  be  represented  by  a  finite  number  of  harmonics;  say  K  e  Z,  2)  The  Nyquist  theorem  is 
satisfied.  Assumption  1  implies  that  the  Fourier  transform  in  (14)  can  be  replaced  by 


Si(eJto)  =  pie\ p(  ]X  / ^(co- kQ)Je\p(-jk(pJJ k 


k--K 


4xf,d 


Vi  =  1,2 


(15) 


\  c  J 

Assumption  2  then  implies  that  -K<ka)o  <K ,  V/:,  -K<k<K  in  order  to  avoid  aliasing  of  the  K 
harmonics.  Since  negative  vibrational  and  rotational  frequencies  have  no  physical  meaning,  it  is 
straightforward  to  note  that  from  the  above  assumptions  coa  e  (0,  k).  Furthermore,  since  the  harmonics 
are  symmetric  about  DC,  and  knowing  a  priori  that  an  MD  target  is  present,  it  is  sufficient  to  identify  the 
first  harmonic  corresponding  to  k=±  1,  i.e.,  \coo\<  71  .  The  criterion  0  <<po  <2  n  is  chosen  for  identify¬ 
ing  <po.  Surprisingly,  the  identifiability  criteria  are  similar,  although  not  identical,  to  the  case  of  sinusoidal 
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parameter  estimation.  However,  the  bounds  on  (Oo  and  <po  are  carrier  independent,  unlike  the  constant 
Doppler  case,  where  the  velocity  bound  was  constrained  by  the  higher  carrier  frequency,  /2. 

The  Bessel  functions  take  arguments  4 Therefore,  any  real  argument,  excluding  zero,  is  valid 
because  Jk( 0)  =  0,  \/k^  0.  Since  negative  values  of  the  maximum  displacement,  d,  have  no  physical 
meaning,  then  d  >0.  Moreover,  in  traditional  sinusoidal  FM,  the  argument  4 x^d/c  is  analogous  to  the 
modulation  index  [11],  which  is  strictly  real  and  positive.  Therefore,  0<d  <oo.  It  is  important  to  note, 
however,  that  it  is  impractical  to  constrain  the  initial  range  of  the  target  to  be  within  the  maximum 
unambiguous  range,  while  allowing  the  maximum  displacement  to  be  infinite.  Therefore,  to  be  consistent, 
we  impose  the  criterion  for  d  as  0 <d  <c/4(f2  -  f{)  =  Ru/2. 

It  is  noted  that  typical  indoor  MD  targets,  such  as  rotating  fans,  pendulums,  etc,  are  essentially  nar¬ 
rowband  FM  (NBFM).  For  example,  consider  the  motion  of  a  human,  the  maximum  displacement  of 
either  the  arm  or  the  leg  is  much  less  than  a  meter,  and  for  a  carrier  in  the  low  GHz  range,  the  bandwidth 
is  typically  a  few  Hz.  In  the  case  of  an  indoor  rotating  fan,  the  received  signal  bandwidth  is  proportional 
to  the  length  of  the  fan’s  blade,  and  is  clearly  a  few  Hz  for  such  carrier  frequencies. 

3)  Accelerating  Target 

The  range  for  this  motion  profile  is  given  by, 

R(n)  =  Ra+  cos (0)(an  +  f3n\  0g  (-/t  /  2,*  /  2)  (16) 

where  the  parameters  a  and  (3  represent  the  initial  velocity  and  half  the  acceleration  (or  deceleration)  of 
the  target,  respectively.  As  before,  R()  is  the  initial  range  of  the  target  and  its  identifiability  is  identical  to 

the  constant  Doppler  case.  Substituting  (16)  in  (4),  we  obtain  the  desired  signal  returns  for  the  accelerat¬ 
ing  target.  It  is  clear  that  an  accelerating  target  represents  a  linear  chirp,  or  a  second  order  Polynomial- 
phase  signal  (PPS)  [20].  The  identifiability  criteria  are  developed  by  putting  0-0  in  (16).  The  model  is 
similar  to  the  chirp  parameter  estimation  model,  albeit  a  few  differences  in  identifiability  [21].  The  terms 
4 Kf^/c  and  4 nf^/c  Vi,  i  =  1,2  represent  the  frequency  and  chirp  rate  in  the  chirp  parameter  estima- 
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tion  model,  respectively.  Hence  for  identifiahility,  along  with  the  fact  that  f2  >/p  we  must  have 
-c/4  f2  <a<c/4f2  and  -c/8  f2  <  /?<  c/8  f2 ,  which  are  carrier  frequency  dependent.  It  is  noted  that  the 
hounds  are  for  the  discretized  versions  of  the  parameters.  For  this  motion  profile,  /zf  (/2;  \|/)  =  0 . 

We  note  that  in  indoor  environments,  the  above  three  motion  profiles  can  occur  individually  or  in 
combination.  While  a  rotating  fan  or  a  clock  pendulum  will  exhibit  purely  MD  motion,  the  motion  profile 
of  a  human  walking  can  be  predominantly  modeled  as  a  combination  of  translational  and  MD  compo¬ 
nents.  In  this  chapter,  however,  we  focus  on  three  “stand-alone”  motion  profiles,  namely,  translation,  MD 
and  acceleration  individually.  Combinations  of  these  motions  are  not  considered. 


2.3  Cramer-Rao  Bounds 

In  this  section,  we  derive  the  CRB  for  the  three  motion  profiles  discussed  previously.  The  parameter  6  is 
assumed  known.  In  Appendix-A,  we  show  that  if  0  is  unknown,  then  the  FIM  becomes  singular.  We 
begin  with  the  generalized  CRB  derivation,  taking  into  account  the  effect  of  the  nuisance  parameters.  The 

parameter  vector  including  the  nuisance  parameters  is  denoted  by  i}  =  [i|/, P\*Pi*G\  tOVS/p+Aw*  where 

\\t  =  [iff\*V/2'V/$-V/f  ]7  •  The  parameters  i/fk,  1  <k<  p  are  the  /^desired  parameters,  some  of  which  are 
present  purely  in  the  phase  of  the  signal  returns.  In  general,  the  range  is  a  function  of  both  n  and  ig  ,  i.e. 
R(n)  -  R{n\\\f).  For  a  complex  Gaussian  vector  x,  with  a  covariance  matrix  C,  defined  in  (7),  the 
Fisher  information  elements  are  expressed  as  [22], 


[FV=Tr 


ac  ac 


a/7,  a/7r 


■  [  3.S  H  ,  3s  || 


+  2Re'^^-  C'1 - , 

3/7,  3/7r| 


,  1  <  (<7,r)  <  p  +  4 


(17) 


where  s  =  [s,  s2]r  = 


f  4tt f  ^ 

(p, +h|(v|/))°exp  j - LR(i|/)  (p2 +h2(vj/))oexp 

c 


:Anfl  . . ^ 


R(m») 


with 


c  )\ 

the  exponential  taken  elementwise,  R(i|/)  and  p(. +  h((i|/)  are  IxN  row  vectors,  whose  /ilh  elements  are 
7?(/i;v|/)  and  pt  +  =  1,2  respectively.  That  is,  the  dependency  on  n  is  suppressed.  The  operator 
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4  o  ’  denotes  the  Hadamard  product  or  element  wise  multiplication.  Following  the  work  of  [23],  and  using 


the  vector  differential  operator,  we  obtain 


ds 

a7 
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3  ,1 
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wxi  J  3  Pi  (p2  +h2(vj/))  J 
where  *•’  denotes  the  generalized  Hadamard  product*,  and  ‘  (  )°_l  ’  denotes  the  Hadamard  division  of  a 
vector  or  a  matrix,  i,e.  if  x  =  [jc,  x2...xN],  x0”1  =  [I  /  x,  ,1  /  x2...  I  /  jcw].  .  In  (17),  the  argument  of  the  trace 
operator  is  zero  except  that  which  corresponds  to  the  Fisher  information  of  the  noise  nuisance  parameters 
namely,  a f,  i  =  l,  2.  Partitioning  the  FIM, 


F  = 


F|,  F.;,l 


1  21 


f22'J 


(19) 


The  desired  CRBs  reside  on  the  diagonal  of  the  inverse  Schur  complement  of  F,  with  respect  to  F22  , 


C/?fi(V)  =  [(F11-F12F2-,F21)-,l 


<n 


(20) 


2  1  ' 


Using  (18)-(20),  it  can  be  readily  shown  that 

F"=2Si 
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\T 


i  V 
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(21) 


F12=F27i=[0p  0,  0p  0,,] 

F22  =  Diag  [2/V  /  a]  INIal  /V/cr4  /V/cr4] 


4  For  a  matrix  A  e  cMxN  =[a[  a2....aA,]  and  vector  xe  CWxl  A  •  x e  cMxN  :=  [a,  °x  a2»x . ° x]  :=  x  •  A 
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where  Diag\  ]  is  a  diagonal  matrix  with  entries  as  given  in  the  brackets  and  the  vector  1^,  is  the  row 
vector  of  size  N  having  entries  all  equal  to  one.  Similarly,  the  vector  0p  is  the  column  vector  of  size  p , 
having  all  zero  entries.  The  CRB  of  the  desired  parameters  is  then  expressed  in  a  compact  form  as 

CRB([y) „)  =  {>,  "'I  (22) 

From  (22),  the  CRBs  for  the  desired  parameters  are  completely  decoupled  from  those  of  the  nuisance 
parameters.  The  decoupling  in  (22)  indicates  that  the  variance  of  the  range  parameters  are  insensitive  to 
the  knowledge,  or  the  information  brought  in  by  the  estimation  of  the  nuisance  parameters,  i.e.  they  are 
unaffected  whether  the  nuisance  parameters  are  known  or  unknown  [22].  From  matrix  Fn  in  (21),  we 
observe  that  the  scalar  multiplier  terms  arise  as  a  result  of  the  carrier  diversity,  implied  by  the  dual 
frequency  operation,  and  therefore  have  an  additive  effect  on  the  net  Fisher  information.  It  yields  lower 
CRBs,  when  compared  to  the  no  diversity  (single  carrier  frequency)  case.  In  fact,  the  scalar  term 

32  n  (/,  /  <J~  +  f2c2)  /  c  indicates  that  the  CRBs  are  the  lowest  when  fx  =  /2  =  0,  in  which  case, 
R  =oo  However,  this  is  a  trivial  case,  as  diversity  ceases  to  exist.  In  general,  since  the  FIM  elements  are 
a  function  of  the  squared  values  of  /,  and  /2,  lower  CRBs  are  achieved  when  higher  carrier  frequencies 

are  chosen.  The  effect  of  the  SNRs  on  the  CRBs  is  straightforward  to  analyze.  The  addition  of  Fisher 
information  is  explored  in  greater  detail  in  the  following  subsections.  In  the  ensuing  analysis,  we  derive 
the  FIM  elements  for  the  desired  parameters. 

2.3.1  Constant  Doppler  CRB 

In  this  case,  the  parameter  vector  is  vj/  =  [R()  vf  .  Substituting  (9)  in  (21)  for  the  constant  Doppler  range 
profile,  we  obtain  the  corresponding  FIM.  In  order  to  make  the  CRB  analysis  more  generic,  we  introduce 
the  parameters,  %:=  f2  /  /,,  and  y/':=a]  /  o\.  It  can  be  readily  shown  that  for  a  single  frequency  radar 
employing  carrier  /, ,  the  FIM,  Fy ,  and  its  inverse  are,  respectively,  given  by  [10] 
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where  K  =  32;r:  is  a  constant.  However,  for  the  proposed  dual  frequency  scheme,  the  FIM  is  given  by 
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It  is  clear  from  (23)  and  (25)  that  the  Fisher  information  for  the  dual  frequency  scheme  can  be  written  as, 

F  =  F/|+F/;  (27) 


The  diagonal  elements  in  (26)  represent  the  CRBs  for  the  initial  range  and  velocity,  respectively.  It  is 
evident  from  the  scalar  matrix  multiplicative  factor  in  (24)  and  (26)  that  the  CRBs  for  both  the  initial 
range  and  velocity  for  the  dual  frequency  case  consistently  assume  smaller  values  as  compared  to  their 
counterparts  in  single  frequency  operation.  Hence,  carrier  frequency  diversity  lowers  the  CRBs  for  the 
range  parameters.  In  (26),  the  matrix  elements  containing  N  characterize  the  temporal  diversity  in  the 
Fisher  information,  whereas  the  scalar  multiplier,  containing  the  SNR  terms  and  carriers,  represents  the 
carrier  diversity  in  the  Fisher  information.  Hence,  the  temporal  diversity  factor  is  identical  for  both 
frequencies,  whereas  the  diversity  induced  by  the  carriers  is  different.  If  X  ~  1,  and  y/' =  1,  with 
pi  ~  p,i  =  1,2,  then  the  carrier  frequencies  are  closely  separated  and  the  SNR  corresponding  to  the  two 
frequencies  are  equal.  In  this  case,  the  net  Fisher  information  is  simply 2F^  and  the  CRB  are  then 
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0.5C7?#(i|/)^  ,  where  CRB(\ |i)^  is  the  CRB  of  the  parameter  vector  vj /,  associated  with  the  single 
frequency  /,. 


2.3.2  Micro-Doppler  CRB 

The  derivation  of  the  CRBs  for  the  MD  motion  profile  is  different  from  that  presented  in  [24],  which  uses 
a  hybrid  combination  of  PPS  and  MD  signals.  The  results  on  the  FIM  elements  presented  in  [24]  are  a 
special  case  of  those  presented  here  and  do  not  provide  closed  form  expressions  of  the  FIM  elements  for 
MD  signals,  neither  do  they  include  the  aspect  angle  and  range  parameters  in  the  employed  radar  model. 
Since  the  frequency  diversity  appends  the  individual  fisher  information  regardless  of  the  motion  profile, 
as  shown  in  (27),  we  can  simplify  notation  by  suppressing  the  terms  involving  the  carrier  frequencies  and 

noise  variances,  i.e. >£u=  2/  af,Vi  =  1,2  and  £2i  —  *  ^  ^  2* 


/ )  MD  Fisher  information 

In  this  case,  the  parameter  vector  is  vj/  =  [Ra  d  (Oa  (po  A p  (pf  The  FIM  can  be  derived  from  the 
general  expression  (21)  using  the  range  equation  (10).  It  is  straightforward  to  show  that 


2  AM 


FR„R„  =  X  X/2 «  (P‘  +  AP‘  COS(COon  -  <t>\  ))2 


1=1  «=0 


2  N- 1 


fr„j  =  cos(<9)^  Yj e2,(P,  +  AP,  cos(a)(ln  -  <pt  ))2  cos(6>„/i  -  <pa ) 


i=  I  0 


2  N -1 


fro0>„  =-cos(d)YJY£2  i(Pi  +&Picos(a)0n-</>i))2dnsm(ct)0n-<p0) 


i=  I  n=0 


2  AM 


frii%  =  cos(0)Y  Y  e2i  ( P .  +  AP,  <20S(C0on  -  </. \  ))2d  sin(6>,n  -  <po ) 


i-  1  0 


FR,AP,  -0'  -0'  FJ4>,  -0 

2  AM 

Fdd  =  ^\6)Y  Y  £2<  (Pi  +  APi cos -  <pt  ))2  cos 2((0on  -  <po ) 


(28) 

(29) 

(30) 

(31) 

(32) 

(33) 
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2  /V-l 


F<t<0o  =~co^(e)Yj'Yj£T->(P>  +  APi  cos(0)on  -  $  ))2  dn  sin(6>0n  -  <po)cos(0)on  -(pa)  (34) 


/=!  n=0 


Z.  /Y  —  I 

Fd%  =cos2(^)^^f2|.(yoi  +  Apicos{(00n-<l>i))2ds\n{(00n-(p0)cos{(00n-(p0)  (35) 
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2  /V-l 
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N-\  N\ 

F*pM=H£ iic°s2(6yj-$),  =^f1(A/72sin2(ft;(,n-^)  (41) 

n=0  n=0 

Since  ht(n\\ \i)  depends  only  on  coo  and  not  on  R0,  d,  and  <po,  the  FIM  elements  in  (32)  are  zero.  In 
fact  the  FIM  elements  for  A px  and  ^  depend  on  0)o  only  and  not  on  the  rest  of  the  parameters.  The  effect 
of  /z-(rc;vj/)  and  hence  A p{  is  explained  in  the  simulations  section.  Further  simplifications  of  the  FIM 
elements  in  eqs.  (28-41)  are  shown  in  Appendix-B.  For  the  special  case  of 
/zy (zi;  v|/)  =  0,=>  Api  =  0,V/  =  1,2  corresponding  to  a  vibrating  target,  the  FIM  is  a  function  of  the  aspect 

angle#,  and  the  desired  parameters  are  v) f~[Ra  d  C0o  <po]T  •  In  order  to  simplify  notation,  the  terms 
with  0  are  separated.  We,  therefore,  obtain  the  final  FIM  as  a  Hadamard  product  of  two  matrices,  the 
matrix  A  containing  the  terms  in  0  and  the  FIM,  F,  corresponding  to  #  =  0. 
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F(0)  =  AoF  = 


(42) 


where 


1  cos(#)  cos(0)  cos(0)  ] 
cos(0)  cos2(0)  cos2(#)  cos2(#) 
cos(0)  cos2(#)  cos2(#)  cos2(0) 
cos(0)  cos2(0)  cos2(#)  cos2(#)J 


(43) 


Using  numerical  inversion  it  can  be  verified  that  F  is  positive  definite  (PD).  However,  from  (43),  A  is 
rank  one  with  one  positive  eigenvalue,  1 +  3cos2(#).  Thus,  A  is  positive  semi-definite  (PSD).  Indeed, 
the  Hadamard  product  of  a  PSD  and  a  PD  matrix  is  a  cause  for  alarm,  and  prompts  us  to  verify  the  non¬ 
singularity  of  F(0),  without  which  the  CRBs  are  unobtainable.  Moreover,  since  F(0)  represents  a  Fisher 
information  matrix,  it  must  at  least  be  PSD,  in  which  case  the  CRBs  are  infinite.  The  positive  definite¬ 
ness  of  F(#)  can  be  proven  using  Lemma-1 ,  implying  non-singularity  of  the  FIM  [25).  In  Lemma-1 ,  the 
matrices  are  assumed  complex,  however  the  results  can  be  generalized  to  real  matrices.  In  [25],  the 
authors  have  provided  the  complete  definiteness  characteristics  of  the  Hadamard  product  of  two  matrices. 
However,  in  our  case,  it  is  sufficient  to  prove  that  F(0)  is  PD. 

Lemma  1 :  Consider  matrices  A,  Be  CVx\  where  A  is  Hermitian  PSD  and  B  is  Hermitian  PD.  Then, 
A  o  B  is  PD  if  all  the  diagonal  elements  of  A  are  not  equal  to  zero. 

Proof :  This  lemma  is  proved  in  [25,  pg.  309,  theorem  5.2.1).  The  proof  uses  the  Sylvester’s  law  of  inertia 
extensively. 


The  CRBs  are  present  along  the  diagonal  of  F(0)  \  From  Lemma-1 ,  F(0)  is  invertible,  we  can  simplify 

nor'  as 
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(44) 


For  deriving  (44),  the  theorem  in  Appendix-C  was  used.  The  inverse  of  the  second  matrix  in  (44)  is 
tedious  and,  therefore,  is  computed  using  numerical  techniques.  It  is  evident  from  (44)  that  the  parameter 

R()  is  insensitive  to  6 ,  whereas  the  other  parameters  are  sensitive  to  0,  depending  on  cos2($). 


2.3.3  Accelerating  Target  CRB 

In  this  case,  the  parameter  vector  is  defined  by  \|>  =  [/?0  a  fi]T .  The  FIM  can  be  derived  from  the 

general  expression  (21)  using  the  range  equation  (16).  Proceeding  in  the  same  fashion  as  before,  i.e., 
using  the  Hadamard  product  approach,  we  obtain 
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The  CRB’s  after  inverting  the  FIM,  are  given  by 


CRB(R0)  = 


3/V2  -3/V  +  2 


{p{£2\  "F  p\& 22)  N(N  +  3/V  +  2) 


C/?fi(y9)  = 


12 


16/V2-30/V  +  ll 

_ X _ 

cos2(0)(p2£2l  +  p\e22)  N(N 4  -5/V2  +4) 


180 


cos2(<9)(/?2£2i  +  p\e22)  N(N4  -  5/V2  +4) 


where  f2i  is  defined  before  and  includes  both  carrier  frequency  and  noise  variance  terms. 


(49) 


(50) 

(51) 


2.4  Simulations  and  Quasi-Experimental  Results 

2.4.1  Numerical  simulations 

We  begin  by  examining  the  effect  of  6  on  the  CRBs.  The  term  l/cos2(#)  in  the  CRBs  will  simply 
increase  the  CRBs  for  increased  value  of  0.  The  CRB’s  become  infinite  when  6  =  ±7y^.  Moreover, 

6  =  ±7y^>  yields  a  singular  FIM.  The  minimum  CRB  is  achieved  when  <9  =  0,  i.e.,  when  the  target  is 

moving  along  the  LOS  of  the  radar  system.  Without  loss  of  generality,  in  the  simulations  we  assume, 
pt  -p,  and  /z,  (/?;  \j/)  =  /i(rt;y),V/  =  1,2.  This  then  implies  Ap(  =  Ap  and  <p{=(p.  The  SNR  for  constant 

Doppler  and  accelerating  target  motion  profiles,  are  denoted  by  SNRi  =  p2  /  a] ,  /  =  1,2  for  the  two 
carriers  respectively.  Fig  1  illustrates  the  sensitivity  of  the  CRBs  for  the  constant  Doppler  and  accelerat¬ 
ing  target  motion  profiles,  the  carrier  frequencies  were  chosen  to  be  /,  =906  MHz  and  f2  =  9 19 MHz. 
Fig.  2  shows  the  CRBs  for  the  MD  motion  profile  for  the  vibration  model  (Ap  =  0),  and  the  rotation 

model.  We  define  SNR  for  MD  as  S7V/?-  =  (p  +  ApY  /  cx  ,  /  =  1, 2 .  In  order  to  be  fair  in  our  comparison, 
the  SNR  for  both  the  vibration  and  rotation  models  was  assumed  identical  in  Fig.  2,  and  the  parameters 
for  vibration  and  rotation  model  are  subscripted  by  V’  and  V’  respectively.  Fig.  2  clearly  shows  that  the 
CRBs  for  the  rotation  model  are  higher  than  those  of  the  vibration  model  for  the  same  SNRs.  It  is  well 
known  that  the  CRBs  of  the  phase  parameters  increase  with  the  bandwidth  of  the  signal  amplitude  [22, 
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26].  The  bandwidth  of  the  amplitude  for  the  vibration  model  is  zero  (DC),  whereas,  bandwidth  of  the 
amplitude  for  the  rotation  model  is  (Oo.  Indeed,  the  scalar  function  h(n;\ |/)  in  the  rotation  model  smears 
the  exponential  part  of  the  returns  and  makes  parameter  estimation  more  difficult.  In  Fig.  1  we  observe 
that  v  has  higher  sensitivity  than  R0 ,  since  it  has  a  lower  CRB  for  increasing  N  [22].  Similarly,  from 

Fig.  1 ,  the  chirp,  or  the  accelerating  target  parameter /?,  is  the  most  sensitive,  whereas  R0  is  the  least 
sensitive.  From  Fig.  2,  for  data  record  length  N  >  50,  the  parameter  coa  is  the  most  sensitive.  It  is 
noteworthy  that  the  CRB  for  Ra  is  much  lower  than  the  CRB  for  parameter  dy  both  of  which  have  the 
same  units  of  measurement  and  represent  a  distance  and  displacement,  respectively.  Sensitivity  of  the 
parameters,  as  evidenced  by  the  CRBs,  are  an  indication  of  some  order  of  estimation  of  the  parameters. 
For  example,  we  have  seen  that  v  is  the  most  sensitive  parameter  for  the  constant  Doppler  motion  and  it 
represents  the  frequency  in  the  sinusoidal  parameter  estimation  problem,  w  hereas  R0  represents  the  phase 
in  the  sinusoidal  parameter  estimation  model.  It  is  well  known  that  the  maximum  likelihood  estimator 
(MLE)  for  the  sinusoidal  parameter  model  estimates  the  frequency  first,  and  then  the  phase  [15].  Similar¬ 
ly,  for  the  PPS  signal  model,  the  highest  order  PPS  coefficient  is  estimated  first,  and  the  lower  orders  are 
estimated  subsequently  [23,  26-28].  In  our  case,  this  agrees  with  (5  being  the  most  sensitive  for  a  second 
order  PPS  signal,  i.e.  the  accelerating  target  returns.  For  the  case  of  MD,  we  have  a  suboptimal  estimator 
already  in  place  in  eqs.  (13-14).  Clearly,  taking  the  Fourier  transform  of  the  MD  returns,  using  the  DC  as 
a  reference,  and  selecting  the  two  nearest  peaks  around  0  Hz  is  the  quickest  way  to  estimate  CO(). 

2.4.2  Quasi-Experimental  results 

In  urban  sensing  applications,  the  rotational  or  the  vibrational  frequency  is  an  important  parameter  for 
differentiating  animate  from  inanimate  targets.  An  indoor  rotating  fan  with  three  plastic  blades  was  used 
as  the  target  for  the  real-data  collection  experiment.  In  order  to  obtain  strong  single  component  returns, 
one  of  the  blades  was  wrapped  in  aluminum  foil.  The  carrier  frequencies  of  the  dual-frequency  radar  were 
chosen  to  be  /,  =906.3  MHz,  and  f2  =919.8  MHz,  and  the  sampling  frequency  was  set  to  1  kHz.  The 
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radar  returns  at  the  two  frequencies  were  simultaneously  collected  for  a  total  duration  of  Is.  The  center  of 
rotation  of  the  fan  was  at  a  distance  of  1.22m  from  the  antenna  feed  point  along  the  LOS.  The  rotation 
speed  of  the  fan  was  measured  to  be  780  rpm,  which  is  in  agreement  with  the  technical  specifications  of 
the  fan,  and  is  considered  to  be  the  true  value  of  the  parameter  cao.  The  objective  is  to  estimate  the 
rotational  frequency  of  the  fan  for  varying  SNRs,  and  compute  its  variance  for  comparison  with  the  CRB. 
The  data  is  at  a  high  SNR  and  is  virtually  noise  free.  In  order  to  vary  the  SNR,  we  add  complex  white 
Gaussian  noise.  Hence  the  term  quasi-expenmental  is  used.  Since  no  estimation  technique  has  been 
advocated  in  this  chapter,  we  use  the  Fourier  domain  analysis,  as  given  by  eqs.  (13-14)  and  extract  the 
frequency  of  the  first  harmonic.  This  procedure  is  definitely  suboptimal,  the  optimal  method  of  maximum 
likelihood  is  not  pursued  here  due  to  its  associated  computational  complexity  [24].  Extracting  the  first 
harmonic  frequency  is  sufficient  to  obtain  the  estimate  of  the  rotational  frequency,  since  typically  indoor 
inanimate  targets  are  NBFM,  see  [24]  for  more  details  on  estimating  harmonics  of  NBFM.  The  FIM  for 
MI)  is  a  function  of  the  parameters  d  and  <po,  which  are  unknown.  Furthermore,  the  amplitudes  of  the 
returns  are  also  unknown.  We  employ  some  suboptimal  schemes  to  estimate  these  and  use  these  in  our 
FIM  calculations.  The  estimation  of  the  unknown  parameters  is  performed  using  the  real  data  without 
adding  artificial  noise,  these  estimates  are  then  considered  as  the  true  parameters.  In  particular,  since  the 
amplitudes  of  the  returns  are  unknown,  and,  in  general,  are  time-varying,  the  amplitudes  can  be  efficiently 
estimated  using  the  second  order  cyclic  moment,  see  [28]  and  references  therein.  On  the  other  hand,  the 
time-frequency  (TF)  distribution  is  particularly  useful  in  extracting  the  parameter  d.  The  maximum  and 
minimum  instantaneous  frequencies  of  the  MD  signal  are  given  by  +2fid0)o  /  c  corresponding  to  the  two 
carrier  frequencies,  respectively.  The  bandwidth  of  the  MD  signal  can  then  be  defined  as  4  ftda)o  /  c.  The 
bandwidth  is  estimated  from  the  original  returns  using  the  spectrogram,  and  since  the  rotational  frequency 
is  known,  computation  of  d  is  straightforward,  this  value  is  then  taken  as  the  true  parameter  value.  The 
return  for  the  first  carrier  was  used  for  calculating  d .  For  estimating  the  phase  ^,,we  use  the  least 
squares  estimator  as  in  [24],  with  K  =  1  using  the  returns  corresponding  to  the  first  carrier  frequency.  For 


71 


this  procedure,  the  returns  were  zero  padded  and  an  8192  point  FFT  was  computed.  These  parameters 
were  then  used  in  computing  the  CRBs.  Fig.  3(a)  shows  the  spectrogram  of  the  original  returns  corres¬ 
ponding  to  the  carrier  /,.  It  also  shows  the  signature  using  the  estimated  d  and  (po  ,  and  the  true  rotational 
frequency  COo  as  a  solid  black  line,  which  closely  follows  the  TF  signature  of  the  returns.  The  maximum 
and  minimum  of  instantaneous  frequency  are  depicted  with  the  dotted  black  lines.  In  Fig.  3(b),  the  raw 
FFT  of  the  returns  corresponding  to  the  carrier  /,  is  shown,  Fig.  3(c)  shows  the  estimated  second  order 
cyclic  moment  for  the  carrier  frequency  /,,  and  Fig.  3(d)  shows  the  real  part  of  the  returns  for  the  first 
carrier.  From  Fig.  3(c)  it  is  clear  that  the  second  order  cyclic  moment  shows  harmonics  at  the  rotational 
frequency  (Oa,  indicating  the  cyclic  RCS  changes,  and  hence  validating  the  cyclic  amplitude  modulation 
we  have  assumed  for  the  MD  rotational  model.  We  are  now  in  a  position  to  compute  the  CRB  and  the 
simulated  mean  square  error  (MSE)  using  Monte  Carlo  simulations.  Fig.  4  shows  the  CRB  for  the 
rotational  frequency  and  the  MSE  of  its  estimates  vs.  SNRs,  100  Monte  Carlo  trials  was  used  in  the 
simulations.  For  each  carrier  frequency,  the  initial  estimate  was  derived  from  a  raw  FFT  search,  and 
subsequently  the  chirp-z  transform  was  used  to  compute  the  final  estimate  of  the  first  harmonic  around 
the  coarse  initial  estimate  from  the  FFT,  the  search  span  of  the  chirp-z  transform  was  8192  points.  This 
estimation  procedure  is  identical  to  that  used  in  [24]  with  K  =  1.  The  final  estimate  of  the  rotational 
frequency  was  the  average  of  these  estimates  for  the  two  carrier  frequencies.  A  strong  and  consistent  bias 
of  the  order  10  '  was  observed  in  the  raw  MSE  for  SNR  >  20  dB,  as  shown  in  Fig.  4,  the  bias-corrected 
MSE  is  also  shown  and  comes  close  to  the  CRB  for  the  rotational  model.  The  strong  bias  maybe  due  to 
the  chirp-z  transform,  necessitating  optimization  methods  to  be  employed  for  estimation,  for  example  the 
non-linear  least  squares.  For  comparison,  the  CRB  for  the  vibration  model  (the  single  frequency  and  dual 
frequency)  are  also  shown,  the  simulated  MSE  is  far  away  from  the  vibration  model  CRBs. 
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2.5  Conclusions 


In  this  chapter,  we  considered  dual  frequency  radar  for  range  estimation  with  application  to  urban 
sensing.  Three  important  and  commonly  encountered  indoor  range  profiles  were  considered,  namely, 
linear  translation,  micro-Doppler  and  accelerating  target  motion  profiles.  The  CRBs  were  derived  for  the 
initial  range  and  key  target  classification  parameters,  namely,  velocity,  acceleration  and  oscillatory 
frequency  of  targets  respectively,  encountering  linear,  accelerating,  and  simple  harmonic  motions.  A 
parametric  model  specific  to  the  dual  frequency  radar  technique  was  developed  for  the  range  profiles, 
incorporating  the  target  aspect  angles.  It  was  shown  that  if  the  aspect  angle  is  unknown,  then  the  FIM 
becomes  singular,  leading  to  unidentifiability  of  the  parameters.  Closed  form  expressions  for  CRBs  were 
derived  for  the  parameters  associated  with  translational  motion  profiles.  However,  for  the  MD  case,  the 
CRBs  were  intractable  in  closed  form,  and  numerical  matrix  inversion  was  pursued.  It  was  also  shown 
that  the  dual-frequency  approach  provides  lower  bounds  as  compared  to  the  single  frequency  counterpart. 
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APPENDIX-A 


We  demonstrate  by  example,  when  0  is  unknown  and  is  a  desired  parameter,  the  FIM  is  singular. 
Specifically,  the  derivation  is  for  the  constant  Doppler  range  profile,  the  derivation  can  be  generalized  to 
any  other  motion  profile.  For  convenience,  and  without  loss  of  generality,  we  assume  =  p.  The 


parameter  vector  sought  after  is  v|/  =  [Ra  v  Q\.  Using  (17)  and  (9),  and  skipping  the  trivial  details 


eN 

gcos(6>)/V(/V-l) 

2 

-ygsin(fl)/V(/V-l) 

2 


F 

rRay 

Fr„91[ 

F  (ff)  = 

Frov 

F 

vv 

Ke 

F 

L  K>e 

F 

rv9 

Fqo  J 

gcos(6>)/V(/V-l)  -v£sin(6Q/V(/V  -1)  1 

2  2 

gcos2(6>)/V(/V-l)(2/V-l)  -v£sin(2fl)/V(/V-l)(2/V-l) 

6  12 

-v£sin(2fl)/V(/V-l)(2/V-l)  gv2sin2(fl)/V(/V-l)(2/V-l) 

12  6  J 


It  can  be  shown  that  the  FIM  is  singular,  implying  non-identifiable  parameters  for  the  chosen  model  [29]. 
This  can,  however,  be  overcome  if  an  additional  carrier  diverse  radar  system,  oriented  differently  from  the 
original  radar  system,  is  used.  This  involves  an  additional  angle  diversity  term  in  the  Fisher  information. 
Moreover,  the  orientation  of  the  two  radars  with  respect  to  each  other  must  be  known  for  identifiability. 
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APPENDIX-B 


In  order  to  bring  the  FIM  elements  for  MD  in  closed  form,  we  need  the  following  which  can  be  derived 
easily  using  elementary  trigonometry,  and  elementary  calculus.  The  series  are  assumed  to  uniformly 
converge  so  that  the  differential  and  summation  operators  can  be  interchanged. 
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We  are  now  in  a  position  to  express  the  FIM  elements  in  closed  form.  The  following  can  be  readily  shown 
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Theorem- 1 :  Let  A,  B  be  two  complex  matrices  e  CA*V.  The  rank  of  A  is  one  with  non-zero  diagonal 
elements  and  the  rank  of  B  is  N.  Then  (A  °  B)"1  =  B  1  °  (A°_l  )T  =  B  1  o  ( \r  )°_l ,  where  A 0-1  =  [a^'  ] . 
Proof:  Since  A  is  rank  one,  it  can  be  rewritten  as  an  outer  product,  given  by 

A  =  xy " ,  A  7  =  y  V ,  with  x,  y  g  CNxl  (C.  1 ) 

The  Hadamard  product  can  then  be  expressed  as 

A  o  B  =  B  o  Xy w  =  DXBD"  (C.2) 
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It  is  also  clear  from  (C.2)  that  AoB  is  non-singular.  At  this  stage  it  is  instructive  to  note  that  if  any 
diagonal  element  of  A  is  zero,  which  is  tantamount  to  any  element  in  (C.3)  being  zero,  then  the  inverse 
of  AoB  does  not  exist.  This  is  primarily  the  reason  for  imposing  the  condition,  that  none  of  the  diagonal 
elements  of  A  are  zero  in  the  statement  of  the  theorem.  The  proof  is  straightforward  from  hereon,  from 
(C.2) 
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Expanding  the  outer  product  and  rewriting  the  elements  in  terms  of  A  ,  we  obtain 
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It  must  be  noted  that  if  A  has  a  rank  greater  than  one,  then  the  inverse  of  the  Hadamard  product  is  quite 
complicated  to  derive.  Fortunately,  such  a  situation  does  not  arise  in  this  chapter. 
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Fig.  2  Ml)  CRBs. 
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Fig.  3  (a)  Spectrogram  of  the  returns,  (b)  FFT  of  raw  returns,  (c)  Second  order  cyclic  moment. 

(d)  Real  part  of  the  returns. 
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WISE 


Fig.  4  Simulation  results  with  respect  to  SNR. 
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3.  Optimal  Waveform  Design  for  Improved  Indoor  Target  Detection  in 
Sensing  Through-the-Wall  Applications 


Abstract 

This  chapter  deals  with  waveform  design  for  improved  detection  and  classification  of  targets  behind  walls 
and  enclosed  structures.  The  target  impulse  response  is  incorporated  in  an  optimum  design  of  the  trans¬ 
mitted  waveform  which  aims  at  maximizing  the  signal  to  interference  and  noise  ratio  (SINR)  at  the 
receiver  output.  The  interference  represents  signal-dependent  clutter  which  along  with  the  wall  degrades 
the  receiver  performance  compared  to  the  free-space  and  zero-clutter  case.  Computer  simulations  show 
sensitivity  of  the  optimum  waveform  to  target  orientation,  but  depict  an  SINR  enhancement  over  chirped 
waveform  radar  emissions  at  all  aspect  angles.  Numerical  electromagnetic  modeling  is  used  to  provide  the 
impulse  response  of  typical  indoor  stationary  targets,  namely  tables,  chairs,  and  humans. 
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3.1  Introduction 


Waveform  diversity  and  design  has  been  recently  investigated  in  many  applications  of  Radar  and  Com¬ 
munications  [  1  ]-[4] .  In  particular,  the  advances  in  adaptive  radar  transmitters  have  permitted  the  use  of 
sophisticated  pulse  shaping  techniques  in  order  to  improve  the  detection  of  point  and  extended  targets  [5], 
[6].  However,  waveform  design  has  not  been  properly  addressed  in  the  context  of  through-the-wall  radar 
imaging  (TWRI)  [7].  TWRI  is  an  emerging  area  of  research  and  development  with  numerous  civilian  and 
military  applications  [37]-[  1 0].  The  ultimate  TWRI  system  objective  is  to  detect  and  classify  animate  and 
inanimate  objects  behind  walls.  This  is  achieved  by  transmitting,  receiving,  and  processing  radio  fre¬ 
quency  signals  over  the  frequency  range  of  [900  MHz-  5  GHz].  Based  on  extensive  experimental  and 
modeling  results,  it  has  been  decided  that  this  frequency  band  provides  desirable  tradeoff  between  wall 
signal  attenuation,  scene  resolutions,  and  antenna  size  [9]. 

TWRI  applications  have  specific  properties  and  challenges  [10]-[15].  Primarily,  the  presence  of 
walls  distorts  the  signal  and  causes  significant  attenuation  and  time  dispersion.  This  dispersion  may 
present  itself  in  the  form  of  reverberations  and  creates  constructive  and  destructive  interference  in  the 
time-domain  which  subsequently  obscures  the  EM  signatures  of  targets  behind  walls.  Further,  an  indoor 
scene,  in  addition  to  humans,  is  populated  by  few  typical  objects,  for  example,  chairs  and  tables  of 
different  possible  sizes  and  shapes. 

Considering  waveform  design  in  TWRI  is  motivated  by  the  fact  that  in  target  detection  and  identifi¬ 
cation,  an  appropriate  choice  of  the  radar  waveform  has  a  considerable  effect  on  the  amount  of 
information  obtained  about  the  target,  especially  in  the  case  of  spatially  extended  targets.  As  such,  the 
design  of  appropriate  waveforms  is  important  for  improving  the  performance  of  TWRI  systems  [16]. 
Indoor  targets  exhibit  different  frequency  responses  when  illuminated  by  a  wideband  signal. 

The  concept  of  matched  illumination  theory  introduced  in  [17],  [21]  considers  designing  the  trans¬ 
mitter-receiver  pair  that  maximizes  the  signal-to-interference  and  noise  ratio  (SINR)  at  the  output  of  the 
receiver  matched  filter  for  optimal  detection  of  targets  with  known  impulse  responses.  It  is  shown  that 
optimum  waveform  design  produces  a  waveform  that  places  most  of  its  energy  into  the  frequency  band  at 


87 


which  the  frequency  response  of  the  target  of  interest  or  the  combined  target-wall  frequency  response 
exhibits  the  highest  peaks,  while  suppressing  the  energy  in  the  frequency  bands  in  which  most  of  the 
clutter  and  noise  power  reside. 

For  target  identification  and  classification,  the  criterion  is  to  design  a  transmitter-receiver  pair  that 
maximizes  the  square  of  the  Mahalanobis  distance  between  the  echoes  from  two  targets  followed  by  a 
bank  of  matched  filters.  The  filter  with  the  highest  output  classifies  the  target. 

In  this  chapter,  we  extend  the  matched  illumination  concept  to  detection  and  classification  of  targets 
behind  walls  and  inside  enclosed  structures.  We  formulate  the  problem  for  generic  targets,  but  focus  on 
human  as  well  as  two  commonly  used  indoor  items,  namely,  a  table  and  a  chair.  We  address  two  different 
cases  for  optimum  waveform  design  in  through  wall  sensing.  In  the  first  case,  we  account  for  the  wall 
transmission  impulse  response,  whereas  in  the  second  case,  the  wall  presence  is  totally  ignored  and  is  not 
incorporated  in  the  design.  The  latter  represents  what  is  referred  to  as  free-space  waveform  design  We 
demonstrate  the  improvement  in  SINR  and  Mahalanobis  distance  when  using  the  optimum  waveform 
over  a  chirp  signal  of  the  same  energy  and  time  duration,  for  both  the  free-space  and  through-the-wall 
cases.  We  also  consider  waveform  design  in  clutter  and  under  range  resolution  requirements.  Waveform 
design  process  with  uncertainty  in  target  orientation  is  also  examined. 

This  chapter  considers  only  a  single  antenna  monostatic  operation.  Consideration  of  multi-antenna 
including  bistatic,  multi-static,  MIMO  operations  is  important  for  the  purpose  of  target  localization  and 
imaging,  but  is  beyond  the  scope  of  this  chapter. 

In  Section  3.2,  we  provide  a  brief  background  on  applying  the  linear-time-invariant  model  to  the  ra¬ 
dar  scattering  problem  as  suggested  in  [16].  We  also  present  an  overview  of  the  matched  illumination  for 
free-space,  and  propose  a  modification  to  matched  illumination,  where  the  wall’s  impulse  response  is 
incorporated  into  the  optimum  waveform  design  process.  Section  3.3  deals  with  the  effects  of  adding  a 
range  resolution  constraint  to  the  optimization  process  and  illustrates  the  trade  off  between  improving  the 
range  resolution  and  enhancing  the  SINR  Section  3.4  addresses  the  issue  of  designing  the  optimal 
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waveform  for  detection  and  classification  under  uncertainty  in  target  orientations.  Various  FDTD  based 
simulations  are  presented  in  Section  3.5,  and  finally  the  conclusions  are  given  in  Section  3.6. 

3.2  Radar  Scattering  Model  and  Target  Impulse  Response 

The  point  target  model  can  reasonably  describe  the  radar  returns  when  employing  narrowband  pulses.  In 
an  idealized  point  target  model,  it  is  assumed  that  the  echo  is  a  scaled  and  delayed  version  of  the  transmit¬ 
ted  waveform.  In  this  case,  the  output  of  the  receiver  filter  matched  to  the  transmitted  waveform 
produced  an  SINR  that  only  depends  on  the  total  energy  of  the  transmitted  waveform.  With  large  signal 
bandwidth,  the  returned  signal  can  be  considered  as  that  scattered  from  several  points  in  an  extended 
region  in  space.  In  this  case,  the  received  echo  is  a  function  of  the  transmitted  waveform  temporal 
characteristics  as  well  as  the  target’s  impulse  response  [16]. 

A  common  way  to  model  an  extended  target  is  to  consider  the  target  as  a  linear  system  and,  thus,  ex¬ 
amine  the  input/output  relationship  of  the  respective  linear  model.  Since  indoor  targets  (except  humans) 
are  usually  stationary,  we  assume  linear  time-invariant  systems  over  the  observation  period.  The  corres¬ 
ponding  system  model  is  shown  in  Fig.  1 . 

If  z(t)  represents  the  electromagnetic  field  illuminating  the  target,  and  s(t)  is  the  scattered  field  at 
an  arbitrarily  point  in  the  far  field,  then 

s(t)  =  z(t)*q(t)  (1) 

where  q{t)  is  the  target  impulse  response  and  (*)  is  the  convolution  operation.  The  same  waveform 
convolved  with  the  clutter  response  W  (/)  yields  the  signal  dependent  clutter  return  c(0  =  z(t)*Wc(t) . 
The  receiver  waveform  is  then, 

r(t)  =  s(t)  +  x(t)  (2) 

where  x(t)  =  c(t)  +  n(t)  represents  the  total  undesired  contribution  of  clutter  and  noise.  The  receiver 
applies  the  matched  filter,  b^^  (t)  ,  to  produce  the  output, 

y«)  =  bmilch(t)*r(t)  (3) 
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The  concept  of  optimum  matched  illumination  waveform  design  was  originally  developed  to  im¬ 
prove  the  detection  of  low  observable  aircrafts  [21].  Here,  this  technique  is  applied  for  target  detection 
and  classification  in  TWRI  applications.  When  the  target  of  interest  is  present,  the  received  signal  consists 
of  the  reflection  of  the  transmitted  waveform  from  the  target,  the  signal  dependent  clutter,  and  additive 
thermal  noise.  In  the  no-target  case,  the  receiver  signal  has  only  a  thermal  noise  component  and  clutter, 
and  exhibits  less  energy  compared  to  the  target  case.  Thus,  in  a  Gaussian  noise  and  clutter  environment, 
in  order  to  improve  target  detection,  we  need  to  enhance  the  receiver  signal-to-interference  and  noise  ratio 
(SINR). 

The  target  response  to  an  UWB  signal  is  obtained  by  full-wave  numerical  simulations  using  the 
commercial  software  package  XFDTD  based  on  the  Finite-Difference  Time-Domain  method  (FDTD). 
Once  the  scattered  field  is  calculated,  an  approximation  of  the  impulse  response  of  a  target  can  be 
obtained  by  deconvolving  the  incident  waveform  from  the  backscattered  field  [22]. 

For  convenience  of  mathematical  representations,  we  use  matrix-vector  formulations.  The  transmit¬ 
ted  signal  z(t )  is  given  by  N  equally  spaced  time  samples  separated  by  an  interval  At ,  i.e., 

Z  =[z0  ,Z|  .  In  a  similar  way,  the  received  signal  s(t)  is  represented  by  M  temporal  samples, 

where  M  >  N  ,  55[j0,j,  ] } .  A  matrix  representation  of  this  model  is, 


^0 

0 

0  ••• 

0 

1 

’  <7o 

1 

■s'fl  "1 

^0 

0  ••• 

0 

si 

*3 

*N-1 

ZN- 2 

ZN 

ZN- 1 

, 

ZM-1 

ZM  -2 

Z-M-N 

J 

4n-  1 

J 

_SM  -1  J 

Z  q  s 


90 


where  Z  is  M  -by-  N  convolution  matrix  containing  the  incident  waveform,  q  is  the  unknown  target’s 
impulse  response,  and  s  is  M  -  dimensional  vector  of  the  scattered  field.  The  target  impulse  response  q 
can  be  obtained  as  the  least  squares  solution, 

q  =  (Z"Z)'z"s  (5) 

where  [U\H  represents  the  Hermitian  operator. 

Once  the  target  impulse  response  is  determined,  its  samples  are  arranged  in  an  MxN  convolution 
matrix  Q  identical  in  structure  to  Z.  The  clutter-free  received  signal  can  be  obtained  by  the  multiplication 
of  the  convolution  matrix  Q  with  the  transmitted  waveform  Z ,  i.e.,  s  =  Qz  .  Accordingly,  the  system 
output,  after  matched  filtering,  is  given  by, 

y = bl,chr = KLhs + b'Lhx  (6) 

where  blt}iaJch  is  the  matched  filter  of  length  M  —  1  .  In  the  sequel,  we  assume  that  the  target’s  impulse 

response  is  known  and  deterministic,  while  the  clutter  Wc  and  the  noise  ny  are  assumed  to  be  independent 
wide  sense  stationary  stochastic  processes. 

The  autocorrelation  matrix  of  the  received  signal  is  Rr  =  E[rr  }  =  Rs+  R  ,  where  £*[•]  is  the  ex¬ 
pected  value  operator  and  Rs  =  P'\ssfl  ]  =  ssff  y  Rx  =  E[xxu  ]  are  the  autocorrelation  matrices  of  the 

signal  and  the  noise  plus  clutter,  respectively.  Here,  Rx  is  the  Hermitian  Toeplitz  matrix  with  its  / -th 
element  given  by 

r,  =  j{Gn(co)  +  Gc.(co)\Z(a))\2}ejl“dco  (7) 

Tt 

where  Gn(CO),  GJCO),  and  |Z(<y)|"  are  the  power  spectral  density  functions  of  the  noise,  clutter,  and  the 
transmitted  waveforms,  respectively.  For  the  case  when  noise  is  stationary  and  white  Gn(CO)  =  N0,  and  in 
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absence  of  clutter,  Gc{cd)  —  0,  matrix  Rr  becomes  diagonal,  with  diagonal  elements  equal  to  the  noise 


power.  The  filter  output  SINR  can  be  expressed  as  [5], 


lW  J?  h 

y _ match _y_matcj\ 

'  hU  /?  h 

match  x  match 


(8) 


The  optimum  matched  receiver  filter  bmaUh  can  be  written  as  the  Weiner-Hopf  equation, 


bmatch  =  R~'s  =  R'Qz  (9) 

Substituting  Error!  Reference  source  not  found,  into  Error!  Reference  source  not  found.,  the  SINR  of 
the  matched  filter  becomes, 

y mu,ch  ~  ™ax  y  =  sH R~ls  (10) 

*match 

Since  the  optimization  is  for  the  purpose  of  obtaining  the  optimal  transmitted  waveform  z,  Eq. 
Error!  Reference  source  not  found,  can  be  rewritten  in  terms  of  the  transmitted  waveform  z  , 

yop,  =  max  y  match  =  ma  \Z”i2z  (11) 

Z  Z 


where  £2  —  Qu RxxQ . 


The  optimum  transmitted  waveform  is  of  finite-duration  and  finite  energy.  When  the  clutter  is 
present,  this  waveform  can  be  obtained  by  using  an  iterative  procedure  described  in  [21  J.  In  the  case  of 
zero  clutter  Gr(CO)  =  0,  as  shown  in  [17],  [21],  the  SINR  y  is  optimized  if  z  is  proportional  to  the 

eigenvector  of  the  autocorrelation  matrix  of  the  target  £2  corresponding  to  the  largest  eigenvalue. 

In  TWRI  applications,  typical  indoor  targets  can  be  a  priori  defined  and  modeled.  Subsequently,  the 
target  optimum  waveform  may  be  designed  off-line  and  be  ready  foremission  with  system  deployment. 

The  model  of  Fig.  I  can  be  extended  to  account  for  the  walls,  as  described  in  Fig.  2.  The  type  of 
wall  significantly  impacts  the  ability  to  detect,  resolve,  and  image  the  targets  behind  it  [9].  In  this 
chapter,  we  apply  computational  electromagnetics  to  find  the  impulse  responses  of  both  homogeneous 
walls  and  simple  targets.  These  responses  are  then  used  to  design  the  optimum  waveforms,  leading  to 
improved  target  detection. 
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The  free-space  matrix  formulation  in  Section  3.2  can  be  modified  to  incorporate  the  wall’s  impulse 
response  into  the  optimization  process.  In  this  case,  the  received  echo  is  expressed  as, 

s(t)  =  z(t)*  w(t)*  q(t)*  w(-t)  (12) 

where  vv’(/)  is  the  wall’s  impulse  response.  Let  u(t)  =  w(t)* q(t)*  w(—t)  be  the  combined  wall/target 
impulse  response.  Note  that  for  symmetric  walls  such  as  homogeneous  ones,  vv(f)  =  vv(— /).  This  com¬ 
bined  response  can  be  arranged  in  a  convolution  matrix  U  having  the  same  form  as  Z  in 
Eq.  Error!  Reference  source  not  found..  In  this  case,  the  combined  wall/target  autocorrelation  matrix 
is/2^  =  UH R~lU .  Accordingly,  and  similar  to  free-space,  the  optimal  waveform  is  proportional  to  the 
eigenvector  with  the  highest  eigenvalue  of  Qw  ,  when  there  is  no  clutter,  otherwise,  the  optimum  wave¬ 
form  can  be  obtained  using  an  iterative  procedure  described  in  [17]. 

In  the  above  discussion,  we  have  assumed  the  wall  characteristics  to  be  known,  or  a  priori  estimated, 
and  we  only  considered  the  wall  transmission  impulse  response  in  designing  the  target  optimum  wave¬ 
forms.  The  estimate  of  the  wall  electric  characteristics  and  EM  signatures  has  been  a  subject  of  many 
recent  publications,  e.g.,  [18],  where  the  wall  parameters,  thickness  and  dielectric  constant,  are  estimated 
and  used  to  model  the  wall  returns,  which  are  then  subtracted  from  the  data.  This  is  an  effective  way  to 
mitigate  the  wall  reflections  and  separate  them  from  those  which  depend  on  the  target.  The  approaches  in 
[19],  [20],  on  the  other  hand,  do  not  estimate  the  wall  parameters,  but  they  rather  remove  wall  reflections 
by  relying  on  similarity  of  wall  returns  measured  from  different  antenna  positions  parallel  to  the  wall. 
Differential  SAR  or  spatial  filtering  methods  can  be  applied  for  this  purpose. 

3.3  Range  Resolution  Considerations 

Optimum  waveform  design  may  be  subject  to  other  constrains,  in  addition  to  maximizing  the  receiver 
SINR.  The  length,  or  time  extent,  of  the  optimum  waveform  can  be  subject  to  imaging  resolution 
requirements.  For  point  targets,  the  range  resolution  is  given  by, 

AR  =  ct/2  (13) 
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where  T  is  the  emitted  pulse  duration  and  c  is  the  speed  of  light.  If  the  extended  target  has  an  impulse 
response  of  L  samples,  then  two  targets  are  resolved  if  they  are  separated  by  M  / 2  samples,  since  the 

convolution  forces  the  echo  to  extend  over  an  interval  of  M  =  L  +  N  —  \  samples,  where  M  >  N  .  The 
parameter  T  in  Kq.  Error!  Reference  source  not  found,  becomes 

Te=j(L  +  N- 1)  (14) 

J  s 

where  fs  is  the  sampling  frequency  of  the  output  signal.  T  herefore,  the  range  resolution  is  determined  by 

both  the  emitted  pulse  duration  and  the  target’s  impulse  response  duration. 

The  presence  of  a  wall  with  an  impulse  response  of  length  K  samples  will  further  degrade  the  range 
resolution,  since  the  two-way  signal  propagation  through  the  wall  must  now  be  accounted  for.  In  this 
case,  the  temporal  interval  in  Eq.  Error!  Reference  source  not  found,  is  given  by, 

r.  =  r,  +  y(K-l)  =  -Jr[(W  +  i-,)  +  2(/f-i)]  (15) 

J  5  J  S 

Since  the  impulse  responses  of  the  wall  and  the  target  depend  on  different  factors,  most  of  which  are  not 
under  the  waveform  designer  or  system  operator  control,  the  only  parameter  that  can  be  varied  to  achieve 
a  certain  range  resolution  is  the  length  of  the  transmitted  waveform  N  .  It  is  clear  from  Eq. 
Error!  Reference  source  not  found,  that  the  minimum  bound  on  the  range  resolution  value,  i.e.,  A Rmin 
is  achieved  when  N  =  1.  In  this  case, 

-(L+2K-2)  (16) 

J s 

For  any  other  value  of  range  resolution,  A R  >  A Rmw  ,  the  corresponding  pulse  width  N  should  at  most 
be, 

N  =  ZL^--L-2(K  +  \)  (17) 

c 
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This  value,  which  represents  the  constrained  length  of  the  optimum  waveform,  should  define  the  number 
of  columns  in  the  matrix  in  eq.  Error!  Reference  source  not  found.. 

3.4  Application  of  Matched  Illumination  Technique  for  Target  Identification 

The  design  of  the  optimized  waveform  that  maximizes  the  probability  of  correct  target  classification  is 
developed  in  [17],  [21]  and  is  similar  to  those  design  techniques  used  for  optimized  detection.  Maximiz¬ 
ing  the  probability  of  correct  classification  is  equivalent  to  the  maximization  of  the  Mahalanobis  distance 
between  the  two  target  echoes, 

n2  ={sa-sp)HR;\sa-sfi)  (is) 

where  sa  and  s ^  are  the  vectors  corresponding  to  the  echoes  from  targets  OC  and  /?,  respectively.  The 

Hermitian  Toeplitz  matrix  Rx  is  the  autocorrelation  matrix  of  the  noise  and  clutter.  Eq.  (8)  can  be 
presented  in  terms  of  the  targets’  impulse  responses, 

n2=z"£2z  (19) 

where  £2  =  (Qa  —  Q ;i)H  Rx\Qa  ~~  Qp)  •  Similar  to  the  no-clutter  case,  the  optimum  discriminating 

waveform  is  proportional  to  the  eigenvector  corresponding  to  the  highest  eigenvalue  of  £2 .  When  the 
clutter  is  present  in  the  scene,  the  optimal  waveform  is  obtained  via  an  iterative  solution  (see  e.g.,  [5], 
[21]).  In  the  presence  of  the  wall,  Qa  and  Qh  are  replaced  by  the  combined  wall/target  convolution 
matrix  presented  in  Eq.  Error!  Reference  source  not  found.. 

3.5  Target  Detection  and  Classification  with  Target  Orientation 
Uncertainties 

The  optimal  waveforms  are  aspect  angle  dependent  since  radar  targets  have  different  impulse  responses 
for  different  orientations.  So,  in  order  to  achieve  optimal  detection  or  classification  for  a  specific  target, 
one  needs  to  know  the  target’s  orientation  with  respect  to  the  radar’s  antenna  a  priori.  This  information  is 
not  always  available  for  the  radar,  especially  when  the  target  is  behind  walls.  One  simple  way  to  over- 
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come  this  problem  is  to  cycle  through  all  possible  waveforms  for  the  given  target  of  interest.  This  may  be 
justified  by  the  stationary  nature  of  the  objects  and  the  feasibility  of  using  a  large  number  of  pulses  over  a 
long  observation  period.  In  this  case,  the  transmitted  waveform  corresponding  to  the  actual  target  orienta¬ 
tion  will  provide  the  highest  matched  filter’s  output  SINR. 

If  the  above  scheme  is  not  possible  due  the  additional  processing  and  delay  in  acquisition  time,  and  if 
only  a  single  waveform  can  be  transmitted  from  the  radar  for  a  given  known  target,  then  two  different 
approaches  based  on  statistical  characterization  or  subspace  decomposition  can  be  pursued  as  an  exten¬ 
sion  or  alternative  to  the  work  presented  in  this  chapter.  In  both  cases,  only  a  single  waveform  is  designed 
for  each  given  target.  In  the  former  approach,  the  waveform  design  is  based  on  the  target  average  radar 
cross-section  (RCS)  behavior,  characterized  by  the  target  frequency  spectrum,  over  angle,  rather  than  its 
frequency  response.  This  approach  can  benefit  from  the  analyses  in  [27,  28].  In  the  subspace  approach, 
eigendecomposition  can  be  applied  to  the  target  impulse  response  matrix  to  reveal  the  dominant  eigenvec¬ 
tors  that  capture  the  main  target  impulse  response  variations  over  angle.  A  combination  of  these 
eigenvectors  can  then  be  used  as  an  equivalent  target’s  impulse  response.  In  the  above  two  approaches, 
we  can  proceed  with  the  waveform  design  based  on  the  average  impulse  response  and  continue  with  the 
same  steps  described  earlier. 

We  have  examined  the  performance  of  designing  a  single  waveform  by  just  averaging  the  target’s 
impulse  responses  over  angles  (orientations).  We  compare,  in  the  Simulation  section,  the  performance  of 
the  averaged  waveform  and  the  optimal  waveform.  Indeed,  the  waveform  obtained  through  the  use  of  the 
new  averaged  target’s  impulse  response  is  suboptimal  and  will  not  provide  the  same  SINR  or  Mahalano- 
bis  distance  improvement  compared  to  the  case  where  the  target’s  orientation  is  known,  or  properly 
estimated,  and  the  waveform  corresponding  to  that  angle  is  transmitted. 

3.6  Numerical  Simulations 

As  a  proof  of  concept,  the  optimum  waveforms  were  constructed  for  the  human  as  well  as  for  the  two 
simple  stationary  targets  that  are  commonly  found  in  an  indoor  environment,  namely,  a  wooden  chair  and 
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a  wooden  table.  The  geometry  and  dimensions  of  the  table  and  the  chair  targets  are  shown  in  Fig.  3.  All 
targets  were  probed  by  a  vertically  polarized  modulated  Gaussian  pulse  (see  Fig.  4).  The  pulse  is  of  1-8 
GHz  extent,  with  almost  uniform  energy  over  the  band  of  1-3  GHz,  which  is  a  common  band  used  for 
TWRI  applications.  It  is  assumed  that  both  transmitter  and  receiver  are  located  in  the  far  zones  of  the 
targets. 

For  the  numerical  modeling,  a  commercial  EM  simulator  XFDTD  by  REMCOM  was  used.  For  each 
target,  we  performed  a  series  of  monostatic  simulations,  computing  the  scattering  field  for  all  azimuth 

angles  (from  0  to  1 80n  with  a  1°  increment)  in  the  horizontal  plane.  The  target  is  located  in  xy  of  the 
coordinate  system,  and  the  zero  degree  angle  corresponded  to  the  X  axis.  The  azimuthal  angle  (p  was 
measured  from  the  positive  x  axis  in  a  counter  clockwise  direction  and  the  elevation  angle  6  from  the 
positive  z  direction  towards  the  xy  plane,  respectively.  In  all  cases  presented  below,  we  used  6  =  90° . 
Thus,  for  each  target,  the  scattered  field  for  180  aspect  angles  was  collected.  The  scattered  fields  for  the 
chair  and  the  table  targets  obtained  at  the  position  of  the  transmitter-receiver  pair  at  cp  —  0°  are  shown  in 
Fig.  5.  The  corresponding  impulse  responses,  obtained  according  to  Section  3.2,  are  depicted  in  Fig.  6. 
Parts  (z)  to  (f)  below  deal  with  the  two  inanimate  objects  whereas  part  (g)  focuses  on  human  target. 

3.6.1  Free-space  illumination 

The  matched  illumination  concept  was  applied  to  a  series  of  impulse  responses,  and  the  optimal  wave¬ 
forms  were  obtained  for  each  aspect  angle  for  the  targets  considered.  A  comparison  between  the 
frequency  response  of  the  target  and  the  optimum  waveform’s  frequency  response  for  the  same  aspect 
angle  is  presented  in  Fig.  7.  This  figure  shows  that  most  of  the  energy  in  the  optimal  waveform  is  concen¬ 
trated  in  one  or  two  narrow  frequency  bands,  corresponding  to  the  target’s  resonant  frequencies  or  to 
those  of  the  highest  frequency  response. 

In  order  to  compare  the  performance  of  different  waveforms,  we  fix  the  input  SNR  to  the  matched 
filter  in  all  cases.  The  noise  variance  may  be  obtained  as, 


97 


2  P 

 max 


jqSAWw/10 


(20) 


where  Pmax  is  the  highest  power  of  the  received  signal  and  SNRjn  —  10 dB .  The  output  SINRs  corres¬ 
ponding  to  the  transmitted  optimum  waveforms  were  compared  to  those  of  a  chirp  waveform  of  the  same 
duration  and  energy.  The  SINR,  computed  for  every  aspect  angle  for  the  chair  target  using  the  optimum 
signal  and  the  chirped  signal,  is  shown  in  Fig.  8.  It  is  evident  that  the  optimum  waveform  provides  a 
significant  improvement  (more  than  10  dB)  over  the  chirp  signal.  The  matched  filter  is  matched  to  the 
expected  echo  rather  than  the  transmitted  waveform.  We  observe  that  the  symmetry  of  the  target  is 
reflected  in  Fig.  8. 


3.6.2  Behind  the  wall  illumination 

The  same  targets  of  part  (a)  are  now  placed  behind  a  homogeneous  concrete  wall  of  8  inches  thick.  The 
permittivity  of  concrete  was  €  -  6s0.  We  find  the  SINRs  at  the  output  of  the  receiver’s  matched  filter  due 

to  a  chirp,  due  to  the  optimum  free-space  waveform,  and  due  to  the  optimum  waveform  for  the  target 
behind  the  wall.  The  results  for  this  part  are  shown  for  the  table  target.  The  optimum  waveform,  in  the 
case  of  the  wall  places  most  of  its  energy  into  the  frequency  bins  where  both  the  wall  and  the  target’s 
frequency  responses  exhibit  peaks  (see  Fig.  9).  It  is  also  clear  from  Fig.  9,  that  the  optimum  waveform 
designed  with  the  wall  accounted  for  performs  better  than  the  two  other  waveforms,  except  at  two  angles, 

namely  60  and  90",  where  its  performance  is  comparable  to  the  free-space  optimal  waveform.  It  is 
interesting  to  note  that  at  these  angles  a  common  peak  exists  between  the  table  frequency  response  and 

the  wall  frequency  response,  which  is  shown  in  Fig.  10  for  (p  =  60° .  In  this  case,  the  optimum  free-space 
waveform,  in  essence,  has  the  same  spectrum  as  the  one  optimized  for  the  wall  case.  At  other  angles,  the 
optimum  waveform  for  the  wall  case  exhibits  a  different  spectrum  than  its  free-space  counterpart.  The 
plots  of  the  table  frequency  response,  the  frequency  response  of  the  optimal  free-space  waveform,  and  the 

frequency  response  of  the  optimal  wall/target  waveform  for  (p  —  30°  are  shown  in  Fig.  11.  It  is  clear 
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from  that  figure  that  the  optimal  wall/target  waveform’s  frequency  response  places  its  energy  into  the 
frequency  that  jointly  resonates  the  wall  and  the  target,  while  the  optimum  free-space  waveform  has  its 
energy  concentrated  at  a  lower  frequency,  where  the  target  only  resonates. 

3.6.3  Target  discrimination 

In  order  to  evaluate  the  performance  of  optimal  waveforms  for  the  discrimination  between  different 
targets  of  interest,  i.e.,  the  table  and  the  chair,  we  examined  the  Mahalanobis  distances  between  the 
target’s  echoes  for  different  transmitted  signals.  We  compare  the  optimal  discrimination  waveforms,  a 
chirp  waveform,  and  the  optimal  detection  waveforms  for  the  chair  and  the  table.  The  result,  depicted  in 
Fig.  12,  clearly  shows  the  advantages  of  using  the  optimal  discrimination  waveform. 

3.6.4  Target  resolution  and  constrained  design 

As  described  in  Section  3.3,  the  minimum  range  resolution  is  a  function  of  the  target’s  impulse  response 
and  the  pulse  duration.  The  target  impulse  response  itself  is  a  function  of  the  dimensions  of  the  target  and 
its  dielectric  properties.  For  proof  of  concept,  we  computed  the  impulse  responses  for  a  metal  chair  and  a 
metal  table  of  the  same  dimensions  as  those  of  the  wooden  targets.  The  impulse  responses  of  the  metal 
chair  and  table  are  presented  in  Fig.  13.  It  is  clear  from  these  figures,  that  the  impulse  responses  of  metal 
targets  are  shorter  as  compared  to  those  of  the  wooden  targets  (compare  Fig.  6  and  Fig.  13).  Therefore, 
the  minimum  range  resolution  for  the  metal  targets  is  also  smaller.  The  minimum  range  resolutions  for 
different  targets  both  in  the  free-space  (K  =  0)  and  behind  the  wall  calculated  from 

Eq.  Error!  Reference  source  not  found.,  are  summarized  in  Table  I. 

A  tradeoff  between  SINR  performance  and  range  resolution  for  a  wooden  table  in  the  free-space  and 
incorporating  the  wall  is  illustrated  in  Fig.  14  and  Fig.  15,  respectively.  One  may  notice  that  as  the  range 
resolution  improves,  the  SINR  corresponding  to  the  optimal  waveforms  decreases.  The  spectra  of  three 

different  optimal  waveforms  for  the  table  at  ^9  =  0°  with  different  range  resolution  constraints  are 
presented  in  Fig.  16.  In  this  case,  the  three  range  resolutions  exceeding  the  minimum  achievable  resolu- 
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tion  (Table  I)  were  used,  namely  2,  3,  and  5  meters.  The  energy  of  optimal  waveforms  is  concentrated 
around  the  resonant  frequency  of  the  target  (Fig.  16).  The  smaller  is  the  range  resolution,  the  broader  is 
the  optimal  waveform’s  frequency  response.  This  leads  to  degradation  in  the  SINR,  since  the  energy  is 
no  longer  confined  within  the  narrowband  around  the  resonance  frequency  of  the  target. 

The  proposed  matched  illumination  waveform  design  technique  can  be  extended  to  handle  targets 
separated  in  cross-range.  This  is  simply  achieved  by  using  spatial  selectivity  and  beamforming  [27]  at  the 
transmit  and/or  at  the  receive,  depending  whether  a  physical  or  synthesized  aperture  is  employed.  In  this 
respect,  for  each  increment  in  the  aspect  angle,  optimum  waveform  based  designed  detector  can  be 
applied.  The  beamforming,  in  addition  to  separating  targets  in  cross-range,  can  also  lower  multipath 
which  are  incident  on  the  sidelobes.  Within  the  same  beam,  constrained  range  resolution  design,  dis¬ 
cussed  above,  can  be  applied  to  enable  resolving  range  multipaths  due  to  the  target.  The  problem  then 
becomes  similar  to  typical  radar  operations  with  two  distinctions  1)  The  transmitted  waveform  is  target 
dependent,  2)  the  range  resolution  is  relatively  relaxed  to  account  for  the  narrowband  property  of  the 
waveform  and  to  utilize  full  prior  detailed  description  of  the  target  under  consideration. 

It  is  noted  that  the  concept  of  matched  illumination  considers  the  detection  of  a  given  target  with 
known  size,  shape,  and  material,  i.e,  impulse  response.  In  this  respect,  once  the  presence  of  the  target  is 
declared  by  the  detector  based  on  the  use  of  the  optimally  designed  waveform,  the  target  information  is 
readily  available  and  no  additional  information  concerning  the  nature  of  the  target  may  be  obtained  via 
high  resolution  imaging. 

3.6.5  Target  orientation  uncertainty 

In  the  case  of  target  uncertainty,  we  design  an  average  waveform,  as  discussed  in  Section  3.5.  It  is  noted 
that  the  Swerling  I  distribution  of  human  RCS  was  already  shown  in  [23].  In  order  to  evaluate  the 
performance  when  using  the  average  waveform,  we  compare  the  waveform  detection  capability  with  that 
of  a  chirp.  The  result  is  shown  in  see  Fig.  17,  which  underscores  the  fact  that  the  average  waveform  is 
“suboptimal”  and,  indeed,  for  some  aspect  angles,  it  looses  performance  against  the  chirp  waveform.  We 
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have  also  implemented  the  stochastic-based  approach  for  waveform  design  which  computes  the  autocor¬ 
relation  matrix  from  the  target  power  spectrum  density  [27].  The  latter  was  obtained  by  averaging  the 
magnitude  square  of  the  frequency  response  of  the  target  with  various  orientations.  The  results  were 
almost  identical,  up  to  numerical  errors,  to  the  average  waveform  shown  in  Fig.  17. 

3.6.6  Effects  of  clutter 

Since  in  TWRI  applications,  we  always  deal  with  a  non-zero  clutter,  we  have  numerically  simulated  the 
presence  of  clutter  for  the  wooden  table  using  the  iterative  algorithm  [21].  The  clutter  and  the  noise  are  of 
independent  samples  and  the  clutter-to-noise  ratio  (CNR)  was  set  to  10  dB.  In  Fig.  18,  the  S1NR  at  the 
output  of  the  receiver  filter  due  to  the  optimum  transmitted  waveform  for  the  clutter,  was  compared  with 
the  SINR  corresponding  to  the  optimum  waveform  for  no-clutter  case  (CNR=0)  and  a  chirp  waveform, 
providing  all  transmitted  waveforms  have  the  same  duration  and  energy.  We  observe  from  Fig.  18  that  the 
optimized  waveform  obtained  via  the  iterative  technique  provides  better  detection  in  the  presence  of 
clutter,  though  the  overall  SINR  improvement  degrades  as  the  CNR  increases.  It  is  noted  that  we  only 
considered  clutter  that  is  modeled  by  equation  (2).  Clutter  arises  from  target-to-target  or  target-to-wall 
interactions  and  multipath  should  be  included  for  a  full  representation  of  all  clutter  sources  in  the  scene. 

It  is  clear  from  Fig.  9,  that  even  with  the  optimum  waveform,  the  table  has  a  different  SNR  in  the 
echo  for  different  aspect  angles.  We  chose  two  aspect  angles,  one  with  a  relatively  low  SNR  (15°  aspect 
angle)  and  another  with  a  higher  SNR  in  the  echo  (90°  aspect  angle).  Figure  19  shows  the  ROC  curves 
corresponding  to  each  of  the  two  aspect  angles  using  three  different  waveforms  at  each  angle.  The  first 
waveform  is  the  optimal  waveform  for  the  table  behind  a  concrete  wall,  whereas  the  second  waveform  is 
the  optimal  waveform  for  the  table  in  free  space,  and  finally  there  is  a  chirp  waveform.  We  can  see  from 
the  ROC  curves  that  the  optimal  waveform  for  the  table  behind  the  wall  has  a  better  probability  of 
detection,  for  a  given  probability  of  false  alarm,  than  the  other  two  suboptimal  waveforms. 
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3.6.7  Human  behind  wall 


We  have  also  employed  XFDTD  for  modeling  a  human  behind  a  cinder  block  wall.  The  human  model 
used  in  FDTD  simulations  is  the  High  Fidelity  Frozen  Male  Body  shown  in  Fig  20.  The  body  consists  of 
5  mm  cubical  FDTD  mesh  cells  and  constitutes  of  23  different  tissue  types,  each  with  its  own  dielectric 
properties.  An  overall  mesh  size  is  136  x  87  x  397  cells.  Although  the  entire  mesh  has  about  4  million 
cells,  it  occupies  only  128  MB  of  RAM  and  highly  optimized  for  numerical  computations.  Cinder  block 
wall  contains  72  cinderblocks  (6  x  12).  Each  cinder  block  is  40.62  x  20.32  x  20.3  cm  made  of  cement 
with  8=7.66  80  and  o=0.06  S/m.  The  XFDTD  cell  size  was  selected  as  4.5mm  as  directed  by  dielectric 
parameter  of  the  object  under  investigation.  The  computational  domain  is  surrounded  by  8  perfectly 
match  layers  (PMLs).  The  padding  between  the  computational  domain  and  the  PML  was  50  spatial  cells 
which  is  more  than  half  wavelengths  away  and  provides  the  sufficient  convergence  of  the  solution. 

The  scene  is  excited  by  a  Gaussian  pulse  with  the  positive  gain  over  the  bandwidth  of  about  0.7- 
3.5GHz.  The  temporal  duration  of  the  pulse  is  251  temporal  steps.  The  scattered  field  is  collected  in 
plane  parallel  to  the  surface  of  the  wall  with  the  stand-off  distance  of  lm.  Such  distance  ensures  that  the 
diffraction  from  the  walls  edges  does  not  have  a  significant  impact  on  the  scattered  field.  The  temporal 
interval  for  the  collection  of  the  scattered  field  was  sufficient  to  have  all  the  secondary  reflections  and 
multiple  path  effect  taken  into  account. 

Table  II  shows  the  SNR  of  the  echo  from  the  human  behind  the  concrete  wall  scenario.  It  can  be 
clearly  seen  from  the  results  that  the  optimal  through  the  wall  waveform  gives  the  highest  SNR,  followed 
by  the  optimal  human  in  the  free  space  waveform  and  finally  the  chirp,  has  the  lowest  SNR. 

3.6.8  Effect  of  wall  errors 

Optimal  through-the-wall  target  detection  or  classification  requires  perfect  knowledge  of  the  wall  trans¬ 
mission  response.  In  practical  situations,  the  wall  characteristics  may  not  be  known  a  priori ,  and  need  to 
be  estimated.  We  evaluate  the  performance  of  the  optimal  detection  waveforms  for  the  wooden  table  in 
terms  of  the  SINR  in  the  presence  of  ambiguities  in  wall  knowledge.  The  same  wall  and  target  parameters 
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of  Fig.  9  are  used.  Figure  21  shows  the  difference  in  SNR  between  the  exact  wall  and  ambiguous  wall 
cases  as  a  function  of  target  aspect  angle  (0-60  degrees)  under  various  error  conditions.  Indeed,  the  SNR 
is  lowered  when  wall  errors  are  present.  For  the  example  presented,  an  average  of  2dB  drop  in  SNR  is 
observed  for  10%  wall  parameter  error. 

3.7  Conclusions 

In  this  chapter,  the  concept  of  matched  illumination  was  considered  for  through  the  wall  radar  imaging 
applications.  An  optimal  waveform  design  technique  was  applied  for  stationary  target  detection  and 
classification.  The  algorithm  was  tested  on  the  two  sample  targets  commonly  found  in  an  indoor  envi¬ 
ronment,  namely,  wooden  table  and  chair,  and  on  a  simulated  human  model  behind  a  concrete  wall.  The 
performance  of  the  optimum  waveform  was  compared  to  a  conventionally  used  chirp  waveform.  The 
result  showed  significant  improvement  when  using  the  optimum  waveform  in  the  SINR  for  both  cases  in 
which  the  targets  of  interest  were  in  free-space  and  behind  the  wall.  The  proposed  technique  was  tested 
in  no-clutter  and  simulated  clutter  scenarios.  In  both  cases,  the  optimized  waveform  outperforms  the 
chirp  waveform.  The  sensitivity  of  the  optimum  waveform  and  the  corresponding  receiver  performance 
due  to  uncertainty  in  target  orientation  was  also  examined. 

The  optimum  waveform  derived  in  this  chapter  was  for  a  single  antenna  monostatic  operation.  The 
work  can  be  extended  to  deal  with  an  imaging  system  with  physical  or  synthesized  aperture,  and  for 
bistatic  and  multi-static  operations  [24-26].  For  these  schemes,  the  variation  of  the  target  RCS  and  its 
correlations  over  different  aspect  angles  and  at  long  standoff  distances  should  be  carefully  examined  prior 
to  waveform  design. 

This  chapter  has  dealt  with  known  targets  and,  as  such,  considered  the  deterministic-based  approach 
for  waveform  design  using  the  concept  of  matched  illumination.  The  corresponding  stochastic-based 
approach  for  waveform  design  based  on  matched  illumination,  discussed  in  [27,  28],  is  a  welcomed  and 
important  extension  of  the  deterministic-based  approach,  in  which  the  target  impulse  response  correlation 
matrix  is  replaced  by  the  autocorrelation  matrix  that  is  derived  from  the  known  target  power  spectrum.  In 
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this  case,  different  target  sizes  within  the  same  class  of  targets,  such  as  chairs,  tables,  or  humans,  give  rise 
to  different  realizations  of  the  impulse  response,  and  can  be  captured  in  a  single  waveform  based  on 
average  statistical  characterization.  Future  work  in  matched  illumination  waveform  design  in  TWRI 
should  consider  the  stochastic-based  approach  for  different  categories  or  classes  of  targets,  and  assess  the 
improvement  in  signal  to  noise  and  clutter  ratio  as  compared  to  conventional  waveforms. 
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Scene  of  observation 


Fig.  1  Block  diagram  of  system. 
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Scene  of  observation 
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Fig.  2  Block  diagram  of  the  system  in  the  presence  of  a  wall. 
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(a)  Wooden  table  (b)  Wooden  chair 

Fig.  3  3-D  models  of  indoor  targets  used  in  simulations.  The  reference  point,  “zero-degree  angle”  was 

chosen  at  the  direction  of  x  axis. 
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Fig.  4  Incident  Gaussian  pulse  (a)  Time-domain,  (b)  Frequency-domain. 
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(a)  Chair 
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Fig.  5  EM  field  scattered  by  the  targets  for  transmitter-receiver  pair  located  at  (p  —  0" 


(a)Chair 


(b)  Table 


Fig.  6  Impulse  responses  for  the  sample  targets  when  transmitter-receiver  pair  is  located  at  (p-  0° 
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Fig.  7  frequency  responses  of  the  optimal  waveform  capturing  the  targets  frequency  band  with  the  highest 

energy  at  (p  =  0° ,  (a)  Chair,  (b)  Table. 
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Fig.  8  SINR  improvement,  for  the  chair  target,  resulting  from  transmitting  the  optimum  waveform  at  each 
aspect  angle  versus  a  chirp  waveform  of  the  same  duration  and  energy  for  all  aspect  angles. 


Fig.  9  SINR  resulting  from  transmitting  the  optimum  through  the  wall  waveform  at  each  aspect  angle  for 
the  table,  the  optimum  free-space  waveform  when  the  target  is  behind  the  wall  and  Finally  a  chirp. 
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Frequency  [GHz]  x  1Q9 

Fig.  10  Spectra  of  optimum  free-space  waveform,  through  the  wall  optimum  waveform,  wall’s  impulse 

response,  and  table’s  impulse  response  at  (p  —  6CT . 


Fig.  11  Free-space  optimal  waveform,  through-the-wall  optimal  waveform,  wall’s  frequency  response, 

and  table’s  frequency  response  at  (p  =  30° . 
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Fig.  12  Mahalanobis  distances  due  to  the  optimum  through-the-wall  discriminant  waveforms,  optimal 
table  detection  waveforms,  optimal  chair  detection  waveforms,  and  a  chirped  waveform. 
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(b) 

Fig.  13  Impulse  responses  for  (p  —  0° ,  (a)  for  metal  table,  (b)  metal  chair. 


1  14 


Fig.  14  Degradation  in  SINK  as  we  vary  the  range  resolution  constraint  of  the  transmitted  waveform,  in 

free  space. 


Fig.  15  Degradation  in  SINR  as  we  vary  the  range  resolution  constraint  of  the  transmitted  waveform, 

when  target  is  behind  the  wall. 
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Frequency  [GHz]  x  1Q9 

Fig.  16  Spectra  of  three  different  optimum  waveforms  with  different  range  resolution  constraint  plotted 

on  top  of  the  table  frequency  response. 
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Fig.  17  SINR  due  to  the  average  waveform  in  case  of  aspect  angle  uncertainty  versus  the  SINR  from 

using  a  chirp. 

Table  I.  Minimum  required  separation  in  range  to  resolve  two  targets. 


Target 

ARmin  in  the  free-space 
(meters) 

ARmin  behind  the  wall 

(meters) 

Wooden  table 

1.86 

2.9 

Wooden  Chair 

1.86 

2.9 

Metal  Table 

0.6 

1.73 

Metal  Chair 

0.6 

1.73 
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Fig.  18  SINR  at  the  output  of  the  receiver  filter  in  the  case  of  the  chair  being  the  clutter  while  the  table  is 

the  target  of  interest. 
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Fig.  19  ROC  curves  when  the  wooden  table  is  behind  a  concrete  wall  and  is  illuminated  by  the  through- 
the-wall  optimum  waveform,  free  space  optimum  waveform  and  the  chirp,  respectively,  and  at  two 

different  target  aspects,  15°  and  90°. 

Table  II.  SNR  of  the  received  signal  from  a  real  human  model  behind  a  concrete  wall,  using  three 

different  waveforms. 


Human  Behind  the  wall 

SNR  using  Optimum  Free 

24.01  dB 

space  waveform 

SNR  using  Optimum  behind 

26.44  dB 

the  wall  waveform 

SNR  using  a  chirp 

20.76  dB 
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Fig.  20  Geometry  of  the  frozen  human  model  used  in  the  simulation. 
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Fig.  21  Difference  in  SNR  between  exact  wall  and  ambiguous  wall  cases. 
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4.  Matched-Illumination  Waveform  Design  for  a  Multistatic  Through-the- 
Wall  Radar  System 


Abstract 

We  present  the  matched  illumination  waveform  design  for  improved  target  detection  in  through-the-wall 
radar  imaging  and  sensing  applications.  We  consider  a  multistatic  radar  system  for  detection  of  stationary 
targets  with  known  impulse  responses  behind  walls.  The  stationary  and  slowly  moving  nature  of  typical 
indoor  targets  relaxes  the  orthogonality  requirement  on  the  waveforms,  thereby  allowing  sequential 
transmissions  from  each  transmitter  with  simultaneous  reception  at  multiple  receivers.  The  generalization 
of  the  matched  illumination  waveform  design  concept  from  a  monostatic  to  a  multistatic  setting  casts  the 
indoor  radar  sensing  problem  in  terms  of  multiple-input  multiple-output  (MIMO)  operations  and  puts  in 
context  the  offering  of  MIMO  to  urban  sensing  and  imaging  of  targets  in  enclosed  structures.  Numerical 
electromagnetic  modeling  is  used  to  provide  the  impulse  response  of  typical  behind-the-wall  stationary 
targets,  namely  tables  and  humans,  for  different  target  orientations  and  at  various  incident  and  reflection 
angles.  Simulation  results  depict  an  improvement  in  the  signal-to-clutter-and-noise-ratio  (SCNR)  at  the 
output  of  the  matched  filter  receiver  for  multistatic  radar  as  compared  to  monostatic  operation. 
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4.1  Introduction 


Design  and  implementation  of  signature  exploitation  systems  is  an  important  area  of  current  research  [1]- 
[6].  The  fundamental  philosophy  behind  such  systems  is  as  follows.  In  typical  radar  operational  environ¬ 
ments,  the  possible  targets  considered  in  the  scene  of  interest  are  known  in  advance.  Sufficient  a  priori 
information  about  the  properties  and  characteristics  of  these  targets,  such  as  the  target  shape,  size,  and 
composition,  is  available.  The  goal  of  the  sensor  system  is  to  detect  the  presence  of  the  target  of  interest 
and  to  determine  its  location.  Signature  exploitation  systems  make  use  of  the  a  priori  information  for 
improved  target  detection  and  classification. 

Utilization  of  target  signature  for  design  of  transmission  waveforms  for  optimum  detection  and  clas¬ 
sification  has  been  recently  investigated  in  many  applications  of  radar,  such  as  air-to-air  and  air-to-ground 
radar  systems  [1],  [4]-[6].  However,  it  has  not  been  utilized  in  the  context  of  urban  sensing  and  through- 
the-wall  radar  imaging  (TWRI).  This  new  emerging  field  of  research  and  development  involves  the 
process  of  remotely  detecting,  classifying,  and  locating  targets  inside  buildings  or  other  structures.  With 
recent  advances  in  both  algorithm  and  component  technologies,  TWRI  is  emerging  as  an  affordable 
sensor  technology  in  both  civil  and  military  settings  [7]-[10].  It  has  been  successfully  sought  out  for 
surveillance  and  reconnaissance  in  urban  environments,  requiring  not  only  the  layout  of  the  building, 
including  types  and  locations  of  walls,  but  also  detection  and  localization  of  both  moving  and  stationary 
targets  within  enclosed  structures  [I  I],  [12].  This  technology  can  also  be  used  in  rescue  missions,  search 
ing  for  fire,  earthquakes,  and  avalanche  victims  and  survivors,  and  behind-the-wall  detection  and 
surveillance  of  suspected  criminals  and  outlaws  [  1 3]-[  1 4].  For  through-the-wall  radar  applications,  in 
addition  to  humans,  there  are  only  a  finite  number  of  objects  that  are  commonly  found  inside  rooms  and 
behind  walls,  for  example,  chairs  and  tables  of  a  few  different  sizes  and  possible  shapes.  As  such,  the 
underlying  application  provides  ideal  premises  for  considering  waveform  design  based  on  signature 
exploitation. 

Target  matched  illumination  theory  is  a  signature  exploitation  concept  that  considers  design  of  the 
optimum  detection  system,  comprising  matched  transmission  waveform  and  ‘matched  filter’  receiver,  for 
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targets  with  known  impulse  responses  [5]-[6],  [15].  More  explicitly,  it  designs  the  transmitter-receiver 
pair  that  maximizes  the  signal-to-clutter-and-noise-ratio  (SCNR)  at  the  output  of  the  receiver  matched 
filter.  Matched  illumination  was  considered  for  the  detection  of  targets  behind  walls  and  in  enclosed 
structures  in  [16]  for  single  antenna  monostatic  operation.  In  this  chapter,  we  develop  matched- 
illumination  waveform  design  for  through-the-wall  applications  using  a  multistatic  radar  system  in  which 
the  system  transmitters  and  receivers  are  spatially  distributed  over  a  wide  area  rather  than  being  housed 
on  the  same  physical  platform.  This  generalization  of  matched  illumination  waveform  design  from  a 
monostatic  to  a  multistatic  setting  puts  into  context  the  offerings  of  multiple-input  multiple-output 
(MIMO)  operations  to  indoor  radar  imaging  and  sensing  [17]-[19]. 

It  is  noted  that  the  process  of  sensing  targets  behind  walls  typically  consists  of  two  stages,  namely, 
the  detection  stage  and  the  classification  stage.  The  latter  is  more  complex  and  requires  processing  off¬ 
line.  This  chapter  treats  the  detection  aspect  of  the  problem.  We  focus  on  widely  separated  transmitters 
and  receivers,  and  exploit  the  variations  in  the  target  impulse  response  as  a  function  of  the  transmitter-to- 
target  aspect  angle  and  the  transmitter-receiver  bistatic  angle  for  improved  through-the-wall  target 
detection.  The  targets  of  interest  are  either  stationary  or  moving  at  very  low  speed,  thereby  allowing 
sequential  use  of  the  transmitters  with  simultaneous  reception  at  multiple  receivers.  We,  therefore, 
develop  the  matched  illumination  detection  system  for  the  multistatic  radar  employing  a  signal  model 
based  on  single  active  transmitters  and  under  the  assumption  of  coherent  processing  of  the  target  returns 
at  the  different  receivers.  The  waveform  for  each  transmitter  is  optimally  designed.  We  analyze  the 
performance  of  the  proposed  multistatic  system  under  both  no-clutter  and  clutter-plus-noise  scenarios  and 
compare  it  to  that  of  single  antenna  monostatic  systems  employing,  respectively,  an  optimal  matched 
illumination  and  a  conventionally  used  chirp  waveform. 

The  chapter  is  organized  as  follows.  In  Section  4.2,  we  develop  the  signal  model  for  the  multistatic 
through-the-wall  radar  system.  The  test  statistic,  the  probabilities  of  detection  and  false  alarm,  and  the 
receiver  operating  characteristic  (ROC)  curves  are  provided  in  Section  4.3.  In  Section  4.4,  we  extend  the 
matched  illumination  waveform  design  technique  for  multistatic  through-the-wall  operation.  Simulations 
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results  based  on  finite-difference  time-domain  (FDTD)  modeling  of  typical  indoor  targets  are  presented  in 
Section  4.5,  and  Conclusions  are  provided  in  Section  4.6. 

4.2  Problem  Statement  and  Signal  Model 

We  consider  the  detection  of  a  single  stationary  extended  target  inside  a  room  with  homogenous  walls 
using  a  through-the-wall  radar  system  that  employs  multiple  spatially  distributed  transmitters  and  receiv¬ 
ers.  The  target  could  be  an  animate  or  inanimate  object.  The  former  includes  a  human,  whereas  the  latter 
encompasses  cache  of  weapons.  The  proposed  system,  unlike  the  monostatic  radar  operations,  exploits  the 
variations  in  the  target  impulse  response  as  a  function  of  the  transmitter-to-target  aspect  angle  and  the 
transmitter-receiver  bistatic  angle  for  improved  performance  [18].  Indoor  targets,  depending  on  their 
orientations,  may  render  monostatic-only  radar  systems  ineffective  in  revealing  the  target  presence. 
Curvatures  and  edges  can  cause  weak  monostatic  target  radar  cross  section  (RCS)  when  viewed  from 
specific  angles. 

For  through-the-wall  application,  the  targets  of  interest  are  either  physically  stationary  or  moving  at 
very  low  speeds  such  that  they  remain  within  the  same  range  gate  and  assume  the  same  impulse  response. 
Therefore,  the  set  of  transmitted  waveforms  need  not  be  orthogonal  to  facilitate  separation  between 
waveforms  arriving  simultaneously  at  each  receiver.  We  leverage  the  target  invariance  property  to 
propose  a  matched  illumination  multistatic  approach  in  which  a  sequential  multiplexing  of  the  transmit¬ 
ters  with  simultaneous  reception  at  multiple  receivers  is  used  [19].  As  such,  a  signal  model  can  be 
developed  based  on  single  active  transmitters. 

4.2.1  Single  Active  Transmitter 

Let  the  active  transmitter  be  located  at  x(=(jtf,yf)  and  the  M  receivers  be  positioned  at 
xrm  -  (x,m  ’  yrm )>  w  =  1 ,  2,...,  A/,  as  shown  in  Fig.  1.  If  z(t)  is  the  energy-limited  transmitted  signal  and 
sm(t)  is  the  target  return  at  the  w-th  receiver,  then 

sJt)  =  qm(t)*z(t)  (l) 
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where  4*’  denotes  convolution  and  qm(t)  is  the  combined  wall/target  impulse  response  corresponding  to 
the  m-th  transmitter-receiver  pair,  and  is  given  by 

qm  (0  =  u(t)  *  %m  it)  *  «(-/)  *  S{t  -  Tjm ).  (2) 

In  the  above  equation,  u(t)  is  the  transmission  impulse  response  of  the  wall,  which  is  assumed  to  be 
known  and  independent  of  the  transmitter  and  receiver  locations.  This  is  a  reasonable  assumption  since,  in 
general,  the  wall  thickness  is  small  compared  to  the  typical  antenna  stand-off  distance,  which  could  be 
10m  or  more  from  the  wall.  In  this  case,  the  differences  in  the  path  lengths  through  the  wall  as  a  function 
of  antenna  position  are  negligible.  The  deterministic  function  £m(t)  is  the  impulse  response  of  the  target 
corresponding  to  the  m-lh  transmitter-receiver  pair,  and  is  presumed  known.  The  discussion  on  how  to 
address  the  case  when  the  target  impulse  responses  corresponding  to  the  various  transmitter-receiver  pairs 
are  not  exactly  known  is  deferred  to  a  later  section.  The  quantity  T]m  in  eq.  (2)  is  the  propagation  delay 
measured  from  the  transmitter  to  the  target  and  then  back  to  the  m-th  receiver.  It  has  contributions  both 
from  the  propagation  delays  through  the  air  and  through  the  wall  [20].  Since  the  target  location  is  not 
known  beforehand,  the  propagation  delay  and  hence  the  combined  wall-target  impulse  response  qm(t) 
should  be  assumed  unknown.  However,  exact  knowledge  of  the  combined  wall-target  impulse  response  is 
required  to  implement  matched  illumination  waveform  design  techniques.  In  order  to  overcome  this  issue, 
range  gating  is  employed  and  the  matched  illumination  waveform  will  be  designed  separately  for  each 
range  gate  within  which  the  target  could  be  present.  That  is,  the  optimum  waveform  depends  not  only  on 
the  target  and  wall  impulse  responses,  but  also  on  the  target  location. 

Let  vv  m(f)  represent  the  clutter  response  for  the  m-th  transmitter-receiver  pair  and  nm{t)  be  additive 
noise.  Both  the  clutter  and  the  noise  are  assumed  to  be  stationary  stochastic  processes,  which  are  indepen¬ 
dent  of  each  other.  The  m-th  received  signal  rm(t)  is  given  by, 

rmit)  =  smit)  +  cmit)  +  nm(t),  Cjt)  =  wcmit)*zit)  (3) 

where  cmit )  is  the  signal  dependent  clutter  return. 
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It  is  noted  that  in  the  above  model,  we  have  only  considered  the  wall  transmission  impulse  response. 
The  radar  return  from  the  wall  itself  arrives  earlier  than  the  return  from  the  target.  The  wall  reflections  are 
assumed  to  have  been  either  resolved  from  that  of  the  target  or  mitigated  using  effective  wall  return 
removal  techniques  such  as  those  recently  proposed  in  [2 1  ]-[23] .  In  [21],  the  wall  parameters  are  esti¬ 
mated  and  used  to  model  the  wall  returns,  which  are  then  subtracted  from  the  target  scene  data.  The 
approach  in  [23],  similar  to  the  work  in  [22],  does  not  estimate  the  wall  parameters  and  instead  prepro¬ 
cesses  the  data  by  applying  a  filter  in  the  spatial  domain  to  notch  the  zero  spatial  frequency  which 
corresponds  to  the  approximately  invariant  signal  returns  from  the  wall  at  different  antenna  locations. 
This  preprocessing  scheme  is  analogous,  in  its  objective,  to  the  moving  target  indicator  (MTI)  clutter 
filter  operation  in  the  time-domain  where  signal  return  invariance  is  over  time,  not  space. 

To  facilitate  simulation  on  a  computer,  we  use  a  discrete-time  formulation  of  the  signal  model.  The 
transmitted  signal  z(t)  is  given  by  Nz  equally  spaced  time  samples  separated  by  an  interval  Ar,  i.e.,  the 

transmitted  signal  vector  is  defined  as  z  =  [z],Z2y">ZN]T .  Target  impulse  responses 
m  =  ,  can  extend  over  different  finite  time  durations.  However,  for  simplicity  of  presen¬ 

tation,  they  are  assumed  to  be  of  the  same  finite  length.  The  combined  wall/target  impulse  response  qm(t) 
is  sampled  at  the  same  rate  as  the  transmitted  signal  to  produce  an  N q  xl  impulse  response  vector 

qm  =  [ q  m i > q  m 2 > •  •  •  * q mN  1  ar|d  the  m_th  received  target  return  sm(t)  is  represented  by  the  return  signal 
vector,  sm=[sl9s29...,sN  ]r,  where  Ns  =  Nz  +  N  - 1 .  A  convenient  matrix  representation  of  the  w-th 
target  return,  free  of  noise  and  clutter,  can  be  obtained  following  a  similar  formulation  in  [  15]  as 

sm  =  Q  J-  (4) 

where  the  N,  x  N,  combined  wall/target  convolution  matrix  Qm  is  given  by 
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Likewise,  by  representing  the  clutter  return  and  the  noise  by  vectors  of  their  respective  temporal  samples, 
the  m-th  received  signal  can  be  expressed  in  vector  form  as 

r  =s  +c  +n  (6) 

where 

r«=[r.i.r.2.-".r«vf,  cm=[cmI,cm2,...,c  ]\  nm=[nmVnm2,...,nmN  ]‘  (7) 


We  collect  the  received  signals  from  the  M  receivers  in  an  MN s  xl  vector  r 

r  =  Qz  +  c  +  n=  s  +  c  +  n 
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(8) 


4.2.2  A  note  on  range  resolution  and  range  gate 

As  mentioned  above,  the  convolution  of  the  transmitted  pulse  of  length  N .  samples  with  the  combined 
wall/target  impulse  response  of  length  N  samples  forces  the  return  to  extend  over  an  interval  of 
N j  =  Nz  +  N  —  1  samples.  For  free-space  propagations,  two  point  targets  can  be  resolved  in  range  if  the 
target  returns  are  separated  by 

T  =  (N,-\)At  (9) 
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where  T  is  the  received  echo  duration  (same  as  the  emitted  pulse  duration).  Accordingly,  for  the  case  of 
extended  targets  behind  the  wall,  the  separation  between  the  target  returns  should  be 

r  =  (Ns  - \)At  =  (Nz  +Nq-  2) At,  ( 1 0) 

In  the  signal  model  of  eq.  (8),  we  assume  that  the  range  gate  is  determined  by  eq.  (10). 

4.2.3  Multiple  transmitters 

Consider  a  multistatic  through-the-wall  radar  system  with  L  transmitters.  When  the  /- th  transmitter  is 
active  and  emits  the  signal  z,,  the  MNS  x  1  received  signal  vector  r,  is  given  by  eq.  (8),  reproduced 
below  with  the  introduction  of  the  subscript  /. 

r,,  =Q,,Z/+C,, +  n,,  =s,,  +c,(  +n /  =  1>2 . L 
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Note  that  we  have  assumed  that  each  of  the  L  transmitted  signals  has  the  same  length  N  and  the  target 
impulse  responses  apparent  to  the  LM  transmitter-receiver  pairs  are  all  of  length  N ’  .  The  values  of  Nz 
and  N  can  be  chosen  based  on  their  possible  maximum  values,  consistent  with  system  and  target 
specifications  and  properties. 

We  combine  the  L  received  signals,  obtained  from  the  sequential  use  of  the  L  transmitters,  to  form  a 
tall  vector  r/wf  of  length  LMNX,  given  by 

r„„  =  Q„„z„„  +  c„„  +  n  =  s„„  +  c„„  +  n 


's 

1 

’s-, 

1 

"n-, 

1 

’c-, 

1 

r„„  = 

s-3 

.  *v  = 

’  Ctot  “ 

>Q,„,  = 

J 

A 

J 

.no 

J 

A 

J 

Q„  o  ...  o  1 

0  Q,  0 


0  0 


Q, 


(12) 


129 


4.3  Detection  Statistic 

In  this  section,  we  formulate  the  through-the-wall  detection  problem  for  the  multistatic  radar  using  the 
proposed  sequential  multiplexing  approach.  We  assume  that  both  the  noise  n{ot  and  the  signal-dependent 
clutter  return  ctot  in  eq.  (12)  are  independent  zero-mean  multivariate  real  Gaussian  processes  with  known 
covariance  matrices, 


k = £KX, } = tf2i lmN,  .  K = E{ clolcl } 


(13) 


where  <J  is  the  noise  variance  and  l[MN  is  an  identity  matrix  of  dimensions  LMNsxLMNs.  More 

explicitly,  the  noise  is  independent  across  the  receivers  from  one  transmission  to  the  next  and  is  identical¬ 
ly  distributed.  Also,  the  signal-dependent  clutter  is  assumed  to  be  uncorrelated  across  the  widely 
separated  transmitters  and  receivers  in  which  case  the  matrix  Rr  takes  the  block  diagonal  form 


Rf  (R^u ,  Rf  1 2  v  •  m  x.w ) 
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(14) 


where  Rc  /  m  =  E{clmc]m}  is  the  Ns  xNs  clutter  covariance  matrix  corresponding  to  the  m-th  receiver  for 
the  1-th  transmission  and  it  depends  on  the  transmitted  signal  zr 

With  respect  to  a  given  range  gate,  we  define  the  null  and  alternative  hypotheses  as 

Hn: 


Clutter  and  noise  only 
Target  plus  clutter  and  noise 


(15) 


The  log-likelihood  ratio  test  is  given  by 

A(rto/)  - |n  y  =* 

A(r,„,)<|nx  => 


H]  is  true 
Hn  is  true’ 


A(r«)  =  ln 


P(rlol/H0) 


(16) 


where  p(r[ol  /  H{))  and  p( rtot  /  //, )  are  the  conditional  probability  density  functions  (pdfs)  of  the  received 
signal  vtot ,  given  the  null  and  alternative  hypothesis,  respectively.  The  parameter  y  is  the  threshold 
which  maximizes  the  probability  of  detection,  while  controlling  the  probability  of  false-alarm.  Given  the 
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signal  model  of  eq.  (12)  and  under  the  assumption  of  Gaussian  distributed  clutter  and  noise,  the  probabili¬ 
ty  density  functions  for  the  null  and  alternative  hypotheses  are 


P(r„„/Ho)  = 

pK,/h  .)  = 


_ [  _ 

{2k)  I  Rr  +<t  ^lmn,  I 

_ I _ 
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exp(  r™(Rr  ■*"  &  I  LMN,  )  rn») 

exp(-(r,„, -s,„,)r(R(  . 
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(17) 


Using  the  expressions  for  p{ rM  /  H0)  and  p(rtol  ///,)  in  eq  (16)  and  absorbing  the  term  independent  of 
r  into  the  detection  threshold,  the  log-likelihood  ratio  test  takes  the  form 


A(r ,„,)>/ 
A(r,J</ 


//,  is  true 
H0  is  true 


,  A  ( r„„ ) =  SL  <  R ,  +  °"2 1  /A»v, )  '  r„„ 


(18) 


with  Y  =  In y+s^,((R(  +a  \LMN  )  'sM.  Thus,  the  optimal  detector  is  a  matched  filter  whose  impulse 


response  bM[1  is  given  by, 

^  match  =  (Rf  ®  ®  LMN,  )  Sln» 

Using  eqs.  (13)  and  (14),  the  matched  filter  of  eq.  (19)  can  be  expressed  as, 
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and  the  output  y<ml  of  the  matched  filter  detector  is  given  by, 

L  M  L  M 

y<>ut  “  ^  S  k  march, l.nflm  “  ®  ^  V,  )  I* /m 

/=!  m=l  /=]  m=l 


(19) 


(20) 


(21) 


It  is  noted  that,  due  to  the  block  diagonal  nature  of  the  matrix  (R  +  <7 lIMN  ),  the  processing  of  the  LM 


received  signals  corresponding  to  the  L  sequential  transmissions  becomes  decoupled.  With  the  1-th 
transmitter  active,  the  corresponding  m-th  received  signal  is  passed  through  a  matched  filter  whose 
impulse  response  entails  the  clutter  statistics  and  the  target  echo  for  the  (1,  m)-th  transmitter-receiver 
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pair  only.  The  LM  matched  filter  outputs  corresponding  to  the  L  sequential  transmissions  are  coherently 
combined  to  obtain  the  output  of  the  overall  multistatic  optimal  detector. 


4.3.1  Receiver  operating  characteristic 

The  probabilities  of  detection  and  false  alarm  for  the  matched  filter  detector  of  eq.  (18)  are  given  by 

Pd  =  J/a(A  /  H})dA,  Pfa  =  J/a(A  /  H0)dA 

y  y 


(22) 


where  /a(A/A/,)  and  /a(A///0)  are  the  likelihood  ratio  distributions  under  the  alternate  and  null 
hypothesis,  respectively.  Since  matched  filtering  is  a  linear  operation  on  the  Gaussian  vector  rtot ,  A  is 
also  Gaussian  distributed  with  respective  means  equal  to  zero  and  bTmatch  sJOJ  under  H0  and  //,,  and 
covariance  b^f(/l(Rc  H-cr2ILW^  )bmaft7r  under  both  hypotheses.  Therefore,  Pd  and  Pfa  take  the  form  [24] 
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\J(bmakh(Rc  +  O  ^IMN,  ^  match  ) 


where  0(w)  =  —  +  —  erf  ~^= 

2  2  \y[2 


zero  mean  and  unit 


is  the  cumulative  distribution  function  of  a  Gaussian  random  variable  with 


it  variance  and  erf  (jc)  =-7=  [e  f  dt  denotes  the  error  function  [25]. 

The  corresponding  receiver  operating  characteristics  curve  is  given  by  [24], 

Pd=(l-<!>(<t>-\\-Pfa)-d)) 

0  hr  5  5 
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b'  s  sT  b 

match  tot  tot  match 


(24) 


(^match(^c  +  ^ LMN  ,  match  ) 


where  d2  can  be  readily  identified  as  the  SCNR  at  the  output  of  the  matched  filter. 
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4.3.2  Coherent  vs.  combining 

Since  the  optimal  detector  requires  coherent  combining  of  the  received  signals  resulting  from  L  sequential 
transmissions,  phase  synchronization  across  the  spatially  separated  transmitters  and  receivers  is  required. 
This  requirement  is  not  too  stringent  for  through-the-wall  applications  in  urban  environments.  For 
illustration,  consider  a  typical  urban  data  acquisition  scenario  shown  in  Fig.  2(a).  Assuming  access  to 
only  the  front  of  a  lOm-by-lOm  building  of  interest,  the  urban  setting  permits  a  maximum  viewing  angle 
of  53°  [26].  For  a  standoff  distance  of  10m,  this  translates  to  a  maximum  distance  of  19.94m  between  the 
sensors. 

However,  in  the  case  when  access  is  available  to  multiple  sides  of  a  building  in  an  urban  setting  or 
for  sensing  in  a  suburban  environment  where  the  buildings  along  a  street  are  widely  separated  (see  Fig. 
2(b)),  the  maximum  distance  between  the  sensors  could  be  relatively  large.  In  such  situations,  it  may  be 
still  feasible  to  distribute  the  transmitters  in  such  a  way  that  only  phase  synchronization  across  the 
receivers  is  required,  enabling  coherent  processing,  while  the  matched  filtered  outputs  corresponding  to 
the  L  sequential  transmissions  are  noncoherently  combined.  In  this  case,  the  test  statistic  is  given  by, 
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y  out. i  * 


v  —  b  r 

Jout.l  match. I  t, 
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where  bma(chl  [^mauh.i.\  ^ match. i. 2 


1 =1 

b match 1  Since  r,  is  Gaussian  distributed  with  means  equal 


to  zero  and  bTmatchJ  s{/ ,  respectively,  under  the  null  and  the  alternative  hypothesis,  and  covariance 
bL*./(Rc,/  +<721  match  J  with  R,./  =  diag  (R(JJ ,  RcJJ RcJm„ )  under  both  hypotheses,  then 

each  yout  l  has  a  central  x  distribution  under  //0,  and  a  noncentral  x2  distribution  under  //, .  However, 
as  the  variance  of  the  pdfs  of  r,  depends  on  the  transmitted  signal  through  s  ,  a  different  pdf  applies  to 


each  youtr  Therefore,  the  derivation  of  the  distributions  of  the  test  statistic  of  eq.  (25)  under  H0  and  H] 

and  the  corresponding  probabilities  of  detection  and  false-alarm  becomes  much  more  involved  and  is 
beyond  the  scope  of  this  work. 
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4.4  Matched  Illumination  Waveform  Design 

The  objective  of  the  matched  illumination  waveform  design  is  to  find  the  transmitted  signal  vector  ztot 


that  maximizes  the  signal-to-clutter-and-noise  ratio  at  the  output  of  the  matched  filter  detector  of  eq.  (18). 
Using  eqs.  (19)  and  (24),  the  waveform  design  problem  can  be  formulated  mathematically  as, 


maxSCNR  =  max 


b  s  s  b 

match  tot  tot  match 


hr  ^ 
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(26) 


Exploiting  the  block  diagonal  nature  of  (R  +cr  )  and  QUJf ,  the  SCNR  simplifies  to 


SCNR  =  £ : z]  ( J  Ql  ( R  +  cr2\ )-'  Q,m  V 

/=1  V  m~l  / 

=  Zz/rQ^Rf./ 


(27) 


Since  the  1-th  term  in  the  summation  on  the  right  hand  side  of  SCNR  in  eq.  (27)  depends  only  on  the  1-th 
transmission,  it  is  clear  that  the  original  waveform  design  problem  of  eq.  (26)  is  equivalent  to  designing 
L  individual  waveforms  z,,z2,...,zL,  such  that 


max z[Q^ (R  ,  +  a2 lMN> y'Q,z„  l  =  \,2,...L 


(28) 


In  the  case  of  zero  or  negligible  clutter,  the  solution  to  the  design  problem  of  eq.  (28)  is  proportional  to 

1  T 

the  eigenvector  corresponding  to  the  largest  eigenvalue  of  Sl{  = — rO^O,  f  1 5],  [16].  On  the  other  hand, 

a  '  1 


when  both  signal-dependent  clutter  and  noise  are  present  and  significant,  the  optimal  waveforms 
zpz2,...,zL  can  be  obtained  iteratively,  as  discussed  in  [15]. 

The  above  analysis  was  carried  out  under  the  assumption  of  known  impulse  responses  of  the  target 
corresponding  to  the  various  transmitter-receiver  pairs.  Note  that  for  a  given  range  gate,  the  target 
impulse  response  at  each  receiver  is  a  function  of  the  incident  angle,  the  reflected  angle,  and  the  transmit- 
ter-to-target  aspect  angle,  which  may  not  always  be  available  when  the  target  is  behind  walls.  One  simple 
way  to  overcome  this  problem  is  to  cycle  through  all  possible  optimal  matched-illumination  waveforms 
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for  the  target  of  interest  given  the  multistatic  sensor  geometry.  The  SCNR  varies  as  a  function  of  the 
target  orientation,  the  incident  and  reflected  angles,  even  when  a  corresponding  optimal  pulse  is  used  for 
each  orientation.  Therefore,  the  transmitted  waveform  corresponding  to  the  actual  target  orientation  and 
angles  of  incidence  and  reflection  will  provide  the  highest  SCNR.  It  is  noted  that  the  symmetry  of  the 
target  can  be  utilized  to  reduce  the  number  of  optimal  waveforms  used  to  interrogate  the  scene. 

4.5  Simulation  Results 

A  multistatic  radar  system,  consisting  of  three  sensors,  was  considered.  Two  different  sensor  functionali¬ 
ties,  listed  in  Table  I,  were  simulated.  In  the  first  case  (Case  I),  sensor  SI  acts  as  a  receiver,  S2  is  a 
transceiver,  and  S3  is  a  transmitter,  whereas  in  Case  II,  the  roles  of  sensors  S2  and  S3  are  reversed.  The 
sensors  are  placed  at  a  standoff  distance  of  10m  from  a  0.2in  thick  solid  concrete  wall  with  a  dielectric 
constant  of  7.66  and  conductivity  of  0.06  S/m,  as  depicted  in  Fig.  3.  A  single  stationary'  target  is  located 
1.3m  behind  the  wall.  For  the  first  set  of  simulations,  a  wooden  table,  shown  in  Fig.  4(a),  was  used  as  the 
target,  with  its  long  side  parallel  to  the  wall.  For  the  second  dataset,  the  model  of  a  standing  human  facing 
the  wall  was  employed.  The  particular  human  model  used,  shown  in  Fig.  4(b),  is  the  “High  Fidelity 
Frozen  Male  Body”  by  Remcom  Inc.  The  human's  phantom  consists  of  23  different  tissue  types,  each 
with  its  own  dielectric  properties.  The  model  is  optimized  for  numerical  computations  and  occupies  only 
128  MB  of  RAM  [27]. 

Since  the  optimal  waveform  design  involves  target  impulse  responses,  free-space  simulations  were 
first  carried  out  using  a  commercial  electromagnetic  (EM)  simulator  XFDTD"  by  Remcom  for  computing 
the  impulse  responses  of  both  the  table  and  the  human  apparent  to  each  transmitter/receiver  pair.  It  was 
assumed  that  both  the  transmitter  and  the  receiver  are  located  in  the  far-field  of  the  target.  Each  target  was 
probed  by  a  vertically  polarized  modulated  Gaussian  pulse,  covering  the  1-8  GHz  frequency  band  with 
almost  uniform  energy  over  1-3  GHz,  as  shown  in  Fig.  5.  Dual-polarized  target  detection  and  imaging 
systems  for  through-the-wall  radar  applications  prove  very  effective,  but  their  consideration  in  the 
underlying  multistatic  waveform  design  problem  has  not  yet  been  investigated  and  it  is  outside  the  scope 
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of  this  chapter.  For  each  transmitter/receiver  bistatic  angle,  the  corresponding  target  impulse  response 
was  obtained  as  the  least  squares  solution  [16], 

S  =  (ZrZr'Zrs  (29) 

where  Z  is  a  convolution  matrix,  identical  in  structure  to  eq.  (5),  containing  the  incident  modulated 
Gaussian  waveform,  and  s  is  the  target  return  vector.  The  transmission  impulse  response  of  the  wall  was 
similarly  computed.  Since  XFDTD*  produced  results  for  a  far-field  scenario  with  plane  wave  excitation, 
the  combined  wall/table  and  wall/human  impulse  responses  were  preprocessed  to  reflect  the  multistatic 
radar  functional  modes  for  Cases  I  and  II. 

Optimal  matched  illumination  waveforms  were  First  constructed  for  the  wooden  table  for  the  zero 
clutter  case.  A  range  gate  of  19.93  ns,  centered  at  the  target,  was  used.  The  target  orientation  was  also 
assumed  to  be  known.  The  SCNR  corresponding  to  the  multistatic  radar  system  using  the  optimal  wave¬ 
forms  are  provided  in  Table  II.  For  comparison,  the  optimal  matched  illumination  waveform 
corresponding  to  monostatic  operation  from  the  transceiver,  was  also  obtained  and  the  corresponding 
SCNR  are  given  in  Table  II.  In  order  to  compare  the  performance  of  the  different  waveforms,  we  used  the 
same  total  transmit  power  and  identical  noise  variance  of  0.001  in  all  cases.  From  Table  II,  we  note  that 
the  multistatic  operation  provides  a  significant  increase  in  the  SCNR  over  the  monostatic  case,  more  so 
for  sensor  functionality  mode  I,  and  subsequently,  a  significant  improvement  in  the  ROC  behavior,  as 
shown  in  Fig.  6.  The  difference  in  improvement  over  the  monostatic  operation  for  the  two  modes  can  be 
attributed  to  the  difference  in  the  monostatic  RCS  of  the  table  corresponding  to  the  aspect  angle  apparent 
to  the  transceivers  in  each  mode. 

In  order  to  gain  insights  into  the  characteristics  of  the  optimum  waveform  relative  to  those  of  the 
wall  and  the  target,  we  plot  in  Fig.  7  the  two-way  wall  frequency  response,  the  table  frequency  responses 
apparent  to  the  two  receivers  with  transmitter  1  active,  and  the  frequency  response  of  the  optimal  wave¬ 
form  transmitted  from  transmitter  1,  all  for  multistatic  sensor  mode  I.  It  is  clear  from  the  figure  that  the 
optimal  waveform  places  its  energy  in  a  narrow  frequency  band  that  jointly  resonates  the  wall  and  the 
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target.  However,  we  note  that  the  energy  of  the  optimal  pulse  is  concentrated  only  where  the  wall  and  the 
target  frequency  response  corresponding  to  receiver  2  exhibit  peaks.  The  reason  for  this  is  two-fold. 
Firstly,  as  clear  from  Fig.  7,  the  peaks  of  the  target  frequency  response  corresponding  to  receiver  I  do  not 
align  very  well  with  the  peaks  of  the  wall  response.  Secondly,  the  impulse  response  of  the  table  corres¬ 
ponding  to  receiver  2  has  a  relatively  small  amplitude  than  that  corresponding  to  receiver  1,  as  shown  in 
Fig.  8. 

Next,  we  constructed  the  optimal  matched  illumination  waveforms  for  the  human  for  the  zero  clutter 
case.  A  range  gate  of  18.5  ns,  centered  at  the  target,  was  used.  The  SCNR  corresponding  to  both  multis¬ 
tatic  and  monostatic  operations  using  the  optimal  waveforms  are  provided  in  Table  III.  Again,  the 
multistatic  system  outperforms  the  monostatic  operation  for  both  sensor  functionality  modes.  The  optimal 
waveform  corresponding  to  transmitter  1  for  mode  I  is  plotted  in  Fig.  9. 

For  both  the  human  and  the  wooden  table  in  the  presence  of  noise  only,  the  output  SCNRs  corres¬ 
ponding  to  the  optimal  transmitted  waveforms  for  monostatic  and  multistatic  operations  were  also 
compared  to  those  of  chirp  waveforms  of  the  same  energy  and  respective  durations.  Similar  to  the 
optimal  waveforms,  the  matched  filter  for  the  chirp  waveform  was  matched  to  the  expected  target  echo,  as 
given  by  eq.  (19),  rather  than  the  transmitted  waveform.  The  SCNRs  corresponding  to  the  chirp  wave¬ 
forms  for  the  table  and  human  are,  respectively,  provided  in  Tables  II  and  III.  It  is  evident  that  the  optimal 
waveforms  significantly  outperform  the  chirp  signal  for  both  monostatic  and  multistatic  operations. 

For  the  case  of  non-zero  clutter,  we  have  designed  the  optimum  multistatic  waveforms  for  the  human 
for  sensor  functionality  mode  I  using  the  iterative  algorithm  presented  in  [15].  The  clutter  and  the  noise 
are  assumed  to  be  of  independent  samples.  The  clutter-to-noise  ratio  (CNR)  was  first  set  to  -10  dB  and 
then  to  lOdB.  Table  IV  provides  the  SCNR  at  the  output  of  the  matched  filter  due  to  three  different 
waveforms,  namely,  the  optimum  transmitted  waveform  for  clutter  and  noise,  the  optimum  waveform  for 
noise  only  case,  and  a  chirp  waveform,  all  of  the  same  duration  and  energy.  For  CNR  of  -10  dB,  we 
observe  that  the  optimized  waveform  obtained  via  the  iterative  technique  provides  better  detection  in  the 
presence  of  clutter,  though  the  overall  SCNR  improvement  compared  to  the  noise  only  case  has  degraded. 
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On  the  other  hand,  when  the  CNR  is  10  dB,  all  transmission  waveforms  yield  similar  SCNR.  This  is  in 
compliance  with  the  degenerate  case  of  significant  clutter,  described  in  [15],  wherein  all  transmission 
waveforms  yield  identical  SCNR.  The  optimal  waveform  corresponding  to  transmitter  1  for  CNR=-10  dB 
is  plotted  in  Fig.  10. 

4.6  Conclusions 

In  this  chapter,  we  presented  a  generalization  of  the  matched  illumination  based  waveform  design  for 
through-the-wall  radar  operation  using  a  multistatic  radar  system  for  improved  detection  of  stationary 
targets.  The  slowly  moving  and  stationary  nature  of  the  indoor  targets  permits  sequential  use  of  the 
transmitters,  rendering  the  use  of  orthogonal  waveforms  unnecessary.  Proof  of  concept  was  provided 
using  EM  modeled  data  for  both  a  human  and  a  typical  indoor  inanimate  object,  namely,  a  wooden  table, 
behind  a  solid  concrete  wall.  The  two  cases  of  noise-only  and  clutter-plus-noise  were  analyzed.  It  was 
assumed  that  the  target  signal  returns  to  a  transmitted  waveform  can  be  coherently  processed  at  the 
different  receivers.  Performance  of  the  optimal  multistatic  waveform  was  compared  to  that  of  the  opti¬ 
mum  waveform  for  single  antenna  monostatic  operation  as  well  as  a  conventionally  used  chirp  waveform. 
It  was  shown  that,  for  the  targets  considered  and  in  the  clutter-free  radar  returns,  the  optimal  multistatic 
waveform  provided  a  superior  performance  and  outperformed  the  others  in  terms  of  the  SCNR  at  the 
output  of  the  matched  filter.  In  the  presence  of  clutter,  the  optimal  multistatic  waveform  provided 
improvements  over  the  other  waveforms.  However,  more  extensive  simulations  need  to  be  performed  in 
order  to  properly  assess  the  performance  of  the  optimal  multistatic  waveform  for  cluttered  scenes. 
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Fig.  1  Through-the-wall  multistatic  sensing  scenario  with  a  single  active  transmitter. 
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Fig.  2  Through-the-wall  sensing  scenarios  in  (a)  urban  environment,  (b)  suburban  environment.  (Adapted 

from  [26]). 


Fig.  3  Layout  of  the  simulated  scene. 
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Table  I.  Functionality  of  the  sensors. 


Operational  Mode 

Transmitter  1 

Transmitter  2 

Receiver  1 

Receiver  2 

Case  I 

S3 

S2 

S2 

SI 

Case  II 

S3 

S2 

S3 

SI 

Fig.  4  (a)  3-D  model  of  a  wooden  table,  (b)  High  Fidelity  Frozen  Human  Model 


Table  II.  SCNR  at  the  output  of  the  matched  filter  for  the  table. 


Target 

Functional 

Waveform 

SCNR  (dB) 

Mode 

Monostatic  Operation 

Multistatic  Operation 

Wooden 

i 

Optimal 

1.2617 

9.6067 

Table 

Chirp 

-7.2687 

-0.4066 

ii 

Optimal 

8.2217 

10.8628 

Chirp 

-3.9222 

0.0437 
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(a) 


(b) 

Fig.  5  (a)  Modulated  Gaussian  Pulse,  (b)  Corresponding  spectrum. 
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Fig.  6  ROC  curves  corresponding  to  the  optimal  waveforms  for  the  wooden  table  for  multistatic  and 
monostatic  operation  under  sensor  functionality  modes  I  and  II. 


Frequency  (GHz) 


Fig.  7  Normalized  Magnitude  Spectra  of  the  two-way  wall  response,  target  impulse  responses  apparent  to 
the  receivers  with  transmitter  1  active,  and  the  optimal  waveform  for  transmitter  1,  for  sensor  operational 

Mode  I. 
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t  (ns) 
(a) 
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Fig.  8  Target  impulse  response  corresponding  to  (a)  Transmitter  1/Receiver  1  (b)  Transmitter  1/Receiver 

2,  for  Mode  I. 


Table  III.  SCNR  at  the  output  of  the  matched  Filter  for  the  human. 


Target 

Functional 

Waveform 

SCNR  (dB) 

Mode 

Monostatic  Operation 

Multistatic  Operation 

Standing 

i 

Optimal 

-4.1907 

1.0279 

Human 

Chirp 

-14.8916 

-10.6437 

ii 

Optimal 

-7.7216 

0.7260 

Chirp 

-14.8094 

-10.6277 
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Fig.  9  Optimal  Waveform  for  Transmitter  1  in  Mode  I  for  no-clutter  case.  The  standing  human  was  used 

as  the  target. 


Table  IV.  Multistatic  SCNR  at  the  output  of  the  matched  filter  for  the  human  in  the  presence  of  clutter 

and  noise. 


Target 

Waveform 

Multistatic  SCNR  (dB) 

CNR=  -10  dB 

CNR=10  dB 

Standing 

Human 

Optimal  -  Clutter  and 

Noise 

-24.9107 

-42.7754 

Optimal  -  Noise  only 

-25.9141 

-43.2472 

Chirp 

-29.6381 

-42.1243 
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Fig.  10  Optimal  Waveform  for  Transmitter  1  in  Mode  I  for  non-zero  clutter  and  noise.  The  standing 

human  was  used  as  the  target. 
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5.  High-Resolution  Through-the-Wall  Radar  Imaging  Using  Beamspace 
MUSIC 


Abstract 

The  MUSIC  algorithm  is  a  high-resolution  direction  finding  technique  which  has  been  successfully 
applied  to  enhance  radar  imaging  in  inverse  synthetic  aperture  radar  (ISAR).  Although  this  signal  sub¬ 
space-based  method  has  proven  effective  when  dealing  with  point  targets  and  high  SNR,  it  may  fail  to 
work  when  directly  applied  to  extended  targets  or  target  returns  of  low  SNR.  The  Beamspace  MUSIC 
(BS-MUSIC),  in  which  the  MUSIC  algorithm  is  applied  to  multiple  beams  is  capable  of  locating  spatially 
extended  targets  in  low  SNR  environments.  We  propose  BS-MUSIC  as  an  image  formation  method  for 
indoor  radar  imaging  problems  and  sensing  through-the-wall  applications.  We  compare  BS-MUSIC 
performance  to  conventional  beamforming  and  element-space  MUSIC.  Imaging  results  with  both  synthe¬ 
sized  and  real  data  demonstrate  the  advantages  of  the  proposed  algorithm  in  depicting  targets  behind 
walls. 
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5.1  Introduction 


High  resolution  radar  imaging  has  recently  gained  much  attention,  as  its  military  and  civilian  application 
areas  are  broadening  to  include  new  challenges  in  urban  sensing  and  homeland  security.  One  application 
that  has  traditionally  made  use  of  high  resolution  radar  imaging  is  the  inverse  synthetic  aperture  radar 
(ISAR)  [1-2].  In  ISAR,  point  targets  are  detected  and  located  by  examining  the  phase  information  induced 
from  changes  in  the  distance  between  the  targets  and  the  antenna.  When  targets  are  approximated  by 
points  in  space,  the  two-dimensional  spatial  spectrum  of  the  scene  embeds  sufficient  information  to 
uniquely  determine  the  targets’  location.  As  such,  high-resolution  spectrum  estimation  methods  such  as 
MUSIC  [3],  Capon  [4],  and  ESPRIT  [5]  can  be  used  as  high-resolution  image  formation  techniques  [2,  6- 
8]. 

Most  recently,  “seeing”  through  obstacles  such  as  walls,  doors,  and  other  visually  opaque  materials, 
using  microwave  signals  has  become  a  powerful  tool  for  a  variety  of  emerging  applications  in  both 
military  and  commercial  paradigms.  Through-the-wall  imaging  (TWI)  has  been  recently  sought  out  in 
behind-the-wall  target  detection,  surveillance  and  reconnaissance  [9].  In  addition  to  surveillance  of 
adversaries,  and  detection  and  monitoring  of  suspected  criminals  and  outlaws,  TWI  technology  is  used  in 
rescue  missions  to  search  for  earthquake  and  avalanche  victims,  and  can  also  aid  fire  fighters  searching 
for  survivors.  In  comparison  to  ISAR  application,  through-the-wall  radar  imaging  faces  the  challenges  of 
counteracting  the  impairing  effects  of  unknown  multiple  wall  characteristics,  and  to  deal  with  small 
enclosed  scenes  with  a  large  variety  of  indoor  targets.  These  targets  may  occupy  several  range  and  cross¬ 
range  cells  or  may  be  confined  to  a  single  resolution  cell.  Similar  to  ISAR,  the  two-dimensional  spatial 
spectrum  of  the  scene  using  collected  TWI  data  can  generated  and  signal  subspace  methods  can  be 
directly  employed  for  resolution  enhancement  [10-11]. 

We  apply  the  MUSIC  algorithm  to  the  outputs  of  spatial  beams,  rather  than  to  raw  sensor  data.  The 
latter  is  referred  to  as  element-space  MUSIC  (ES-MUSIC).  Although  different  in  some  ways  from 
conventional  one-dimensional  beamspace  MUSIC,  we  will  refer  to  the  proposed  method  as  beamspace 
MUSIC  (BS-MUSIC),  so  as  to  distinguish  it  from  the  ES-MUSIC.  In  the  underlying  indoor  radar  imaging 
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application,  beamspace  processing  has  three  main  advantages  over  element-space  processing.  First,  BS- 
MUSIC  outperforms  ES-MUSIC  in  low  signal-to-noise  ratio  (SNR)  cases  due  to  the  processing  gain 
offered  by  beamforming.  The  beamforming  gain  is  equal  to  the  number  of  sensors  when  a  single  narrow- 
band  signal  is  used.  This  gain  is  important  in  view  of  significant  signal  power  attenuations  that  can  be 
caused  by  exterior  as  well  as  interior  walls.  Second,  the  ES-MUSIC  relies  on  the  basic  assumption  of 
point  targets  [1].  Extended  targets  violate  the  point-target  assumption  and  cannot  be  properly  processed, 
and  furthermore,  their  presence  prevents  the  ES-MUSIC  from  locating  even  point  targets.  With  no 
assumptions  on  target  characteristics,  in  terms  of  type  and  size,  the  BS-MUSIC  is  considered  attractive 
since  it  can  image  both  point  targets  and  extended  targets  simultaneously.  This  capability  is  essential  for 
indoor  imaging  where,  due  to  small  standoff  distance  and  limited  bandwidth  and  aperture,  many  targets 
behind  walls  classify  as  spatially  extended  targets.  Third,  the  ES-MUSIC  requires  two  dimensional 
interpolation  to  provide  sampled  data  along  a  rectangular  grid.  In  BS-MUSIC,  the  data  sampled  along  a 
rectangular  grid  is  obtained  by  beamforming,  rather  than  interpolation.  Avoiding  interpolation  is  a 
welcomed  step,  given  the  fact  that  a  single  outlier  in  the  data  set  can  even  result  in  large  interpolation 
errors.  It  is  noted  that  rectangular  grid  sampling  in  both  algorithms  is  required  to  permit  covariance  matrix 
estimation  from  one  set  of  measurements. 

We  examine  in  this  write-up  the  BS-MUSIC  performance  with  real  data  as  well  as  synthesized  data. 
The  real  data  is  collected  by  the  series  of  experiments  carried  out  in  the  Radar  Imaging  Lab  at  the  Center 
for  Advanced  Communications,  Villanova  University.  The  test  results  show  that  the  BS-MUSIC  provides 
better  imaging  quality  compared  to  the  ES-MUSIC  and  conventional  beamforming. 

5.2  High-Resolution  Imaging 

5.2.1  Signal  Modeling 

A  narrowband  signal  with  frequency  /()  transmitted  and  received  at  an  antenna  located  at  (xa,ya)  can  be 
modeled  as. 


151 


(1) 


p- 1 

z(/0)  =  Y?p  1  c } 

n= o 

where  a p  is  the  reflection  coefficient  of  the  p  -th  target,  or  the  target  cross-section,  c  is  the  propagation 
speed,  P  is  the  number  of  targets  in  the  scene,  and  r  is  the  range  from  the  antenna  to  the  p  -th  target 
such  that 

rP  =  .  (2) 

where  (x p,yp)  is  the  location  of  the  p  -th  target.  All  targets  assume  fixed  positions  over  the  imaging 

interval.  This  interval  includes  data  acquisition  time  over  all  frequencies,  and  in  the  case  of  synthetic 
aperture  system,  accounts  for  the  time  to  move  the  antenna  from  one  position  to  another.  In  equation  (2), 
the  range  2 rp  implies  that  both  transmit  antenna  and  receive  antenna  are  located  at  the  same  position 

(monostatic  radar).  In  order  to  obtain  a  high  range  resolution  image,  a  wideband  signal  must  be  used. 
Suppose  that  S(f)  is  the  Fourier  transform  of  the  wideband  signal.  Then, 

s(t)=  \  S(f)exp{j2xft}df.  (3) 

A  stepped-frequency  signal  is  an  approximation  of  a  wideband  signal  with  finite  number  of  narrowband 
signals  whose  center  frequencies  are  uniformly  spaced  within  the  signal  bandwidth.  The  stepped- 
frequency  implementation  of  high  resolution  imaging  is  considered  attractive  due  to  its  hardware  sim¬ 
plicity,  function  reliability,  and  ability  to  control  the  signal  power  at  different  frequencies  for 
counteracting  wall  signal  dispersion  [12].  Using  Nf  frequencies,  the  stepped-frequency  signal  is  given 

by, 

N/ 

s(t)  =  '£jS(fm)exp{j27rfmt},  (4) 

w=0 

with 
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(5) 


fm 

A/ 


fL+mAf, 

Jh  ~  jo 

Nf  “1  ’ 


where  the  frequency  band  is  bounded  by  fL  and  fH  .  In  a  monostatic  radar  system  with  Na  antennas,  the 
received  spatio-temporal  signal  assumes  the  two-dimensional  form, 

p- 1 

z(m,n)  =  ]>y( exp[-j4nfmrpit  / c}.  (6) 


In  the  above  equation,  the  frequency  index  is  m  =  -1  and  the  spatial  index,  rp  n ,  is  the  distance 

between  the  p- th  target  and  the  /2-th  antenna  for  n  =  0,...,Afa  -1 . 

For  long  standoff  distances  where  the  center  of  the  scene  is  far  from  each  antenna  in  the  imaging  sys¬ 
tem,  rpn  can  be  approximated  by, 


rP.n  =rn+xP  cos  °n  ~  V,  sin  0„  (7) 

where  rn  is  the  range  from  the  n  -th  antenna  to  the  center  of  the  scene  (x c>yc),  (xpyyp)  is  the  position  of 

the  p  -th  target  in  jc  —  y  plane.  The  center  of  the  scene  is  at  (0,0) ,  and  0n  is  the  looking  angle  of  the  n  - 
th  antenna  (see  Fig.  1  ).  Accordingly,  the  range  r  assumes  a  fixed  value,  independent  of  target’s  location. 
We  note  that  rn  can  be  easily  determined  if  the  scene  represents  an  interior  of  an  enclosed  structure,  such 
as  a  room  or  a  building,  whose  dimensions  are  known  from  using  ground-based  monitoring  or  surveil¬ 
lance  system.  Aerial  mapping,  such  as  UAV,  can  also  provide  accurate  2-D  extent  of  the  structure  and  its 
center.  If  the  phase  term  due  to  r  is  compensated  for  in  each  antenna,  then  equation  (6)  becomes, 

p  i 

z(m, n)  =  exp{-7'47r/m  (xp  cos  On  -  yp  sin  0n )  /  c} .  (8) 

p=0 

The  above  expression  is  exactly  the  same  as  the  1SAR  signal  model  [2].  Let  Z(kx,kx)  be  the  2-D  spatial 
spectrum, 

p- 1  {  Att  .1 

z(Mv)= Z  expi -J — (k*xP  -  k*yP )  r  • 

p= 0  l  c  ) 
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Then,  equation  (8)  can  be  viewed  as  a  sampled  version  of  equation  (9),  i.e., 


z(m,n)  =  Z(kx, 


^  '  1=  fcosO,k=fs\nd  ’ 

x  m  n7  \  *  m  n 


(10) 


v 


Fig.  1  Geometry  of  the  scene  and  the  antennas.  The  x-y  plane  represents  global  locations  whereas  the 

x  -  y  plane  represents  a  scene  with  center  at  (0,0) . 


5.2.2  MUSIC  Based  Imaging 

Estimating  the  covariance  matrix  is  necessary  for  performing  the  MUSIC  algorithm.  The  eigen  decompo¬ 
sition  of  the  covariance  matrix  yields  both  the  signal  subspace  and  the  noise  subspace.  The  latter  is  used 
to  define  the  projection  operation  for  constructing  the  MUSIC  spectrum,  in  which  peak  positions 
represent  an  estimate  of  target  locations.  In  temporal  processing,  estimating  the  covariance  matrix  is 
achieved  by  time-averaging  of  the  instantaneous  correlation  values.  The  main  challenge  in  radar  image 
applications  is  that  correlation  averaging  is  not  applicable,  as  only  a  single  snapshot  (data  set)  is  available. 
However,  it  is  possible  to  estimate  the  covariance  matrix  from  virtual  snapshots  generated  from  the  data 
set  [2,  6].  In  order  to  obtain  the  virtual  snapshots,  the  data  in  (10)  should  be  uniformly  sampled,  which 
can  be  achieved  through  2-D  interpolation.  Linear,  polynomial,  or  spline  interpolation  methods  can  be 
applied.  For  simplicity,  we  used  linear  interpolation  in  the  simulations  in  Section  5.4. 

Let  z(m,n)  be  the  uniformly  sampled  version  of  (10)  for  m  =  0,...,A/  —  1  and  n  =0,...,/V -1  .  That 
is. 


154 


(11) 


z(m,n)  =  £cr  expS -j—(kmxp  - kny'p) 
/>=()  l  C 


where  the  sampled  spatial  frequencies  are 


l  =  k,.»+rnAkx, 


kn  =  ky0  +  nAkf. 


(12) 

(13) 


The  frequency  steps  A kx  and  A ky  in  equations  (12)  and  (13)  represent  the  interpolation  intervals.  Note 
that  the  spatial  frequencies  km  and  kn  are  hounded  by 


(14) 


fo  sin  9a  <  kn  <  f0  sin  6  . 

a 

In  order  to  avoid  aliasing,  the  following  inequality  should  hold. 

-1/2<(A/  -\)Akx,(N -\)Aky  <1/2.  (15) 

Define  Z  as  an  MxN  matrix  whose  m/i -th  element  is  z(m,n) .  Let  F  be  the  KxL  sub-matrix  of  Z 
such  that 

Zj,j  *  S,y+L-1 

p  ^i+\,j  ^i+1  ,j+\  S+l,;+H 

Kj=  :  :  : 

_Zi+K-\,j  Zi+K-\,j+\  Zi+K  1,y>L-lJ 

where  z,  is  the  ij  th  element  of  the  matrix  Z  .  Then,  the  virtual  snapshots  z.  is  defined  as 

Z-  =  V£c{F/(  } 

where 

vec{[\x  •••  vv]}  =  [v"  •••  v"]". 

Let  u ,(xp,yp)  be  the  i  -th  virtual  snapshots  when  a  point  target  is  located  at  (xp,yp).  Then, 


(16) 


(17) 
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p-\ 

z;  =  W,.y,) 

p=i) 

P-\ 

=  +Mv.vP)}u«(^,->V) 

P= 0 

which  shows  only  a  constant  phase  difference  for  each  target  between  consecutive  snapshots. 
The  covariance  matrix  estimate  is  provided  by  averaging  the  outer  products, 

i  w 

r=—  y  z,.zf, 

N,&“ 


(18) 


(19) 


where  Ns  =  m\n[M  -  K  +  1, vV-L  +  1 }  is  the  number  of  virtual  snapshots.  Processing  of  the  virtual 

snapshots  according  to  (19)  amounts  to  subarray  and  subband  averaging,  which  utilizes  the  phase  differ¬ 
ence  to  restore  the  rank  of  the  covariance  matrix  in  coherent  environments  [6].  If  forward-backward 
smoothing  is  used,  the  covariance  matrix  estimate  becomes  [6] 


R  = 


Z  <z-z.w +JzXJ). 


(20) 


where  (  )*  and  (  )r  denote  complex  conjugate  and  transpose,  respectively,  and  J  is  the  exchange  matrix, 
such  as 


0  •••  0  l] 

0-10 

1  0  •••  0J 


Rather  than  using  diagonal  sub-matrices,  as  implied  by  equation  (20),  it  is  possible  to  use  all  sub-matrices 
of  Z  [13].  In  this  case, 


R  = 


1 


M-KN-L 


2(M  -K  +  \)(N-L  +  \)  UU 


Z  Z<ZMZ"  +KX;-D 


(21) 


where  z. :  .  =vec{F  i).  Although  there  are  some  cases  where  the  covariance  matrix  estimate  in  (21)  is 
invalid,  those  cases  are  unlikely  to  occur  in  the  applications  where  the  number  of  targets  P  are  typically 


much  smaller  than  M  and  N  .  The  sufficient  condition  for  the  covariance  matrix  estimate,  R  to  be 
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valid  is  L>  px  and  K  >  py  where  px  and  pv  are  the  maximum  multiplicity  of  the  target  position  along 


x-axis  and  y-axis,  respectively  [14].  Multiplicity  denotes  the  number  of  targets  that  have  the  same 
position  along  one  axis.  For  example,  if  there  are  three  targets  that  have  the  same  x  position,  the  multip¬ 
licity  along  x-axis  becomes  three. 

The  MUSIC  algorithm  based  on  the  covariance  matrix  estimate  in  (21)  has,  in  general,  better  per¬ 
formance  than  that  associated  with  the  estimate  in  (20).  This  is  because  the  errors  in  the  perturbation 
errors  of  the  matrix  elements  are  inversely  proportional  to  the  number  of  samples  [15].  The  latter  is  much 
larger  in  (21)  than  in  (20).  For  example,  for  M  =  N  =  2K  =  2 L,  the  number  of  snapshots  associated  with 
equation  (20)  is  2(^  +  1)  whereas  the  corresponding  number  in  (21)  is  2(K  +  \)2. 

The  MUSIC  estimator  of  target  locations  relies  on  the  eigendecomposition  of  the  covariance  matrix 
estimate  and  is  given  by, 


,,  ,  y)a(x,y) 

aH(x,y)W’\\"a(x,y) 

where  the  virtual  array  response  vector  a(jc,  y)  is  of  length  KL  and  given  by. 


(22) 


a(jc,y)  =  [l  a  •••  aK  1  p  afi  •••  aK  {/3 


J32  a}p2  aK  ]fiL  ']  (23) 

4/r 

a  -  exp  {-j  —  Akxx] 
c 

P  =  exp{;'— Mvy) 
c 

The  matrix  W  in  equation  (22)  is  the  noise  subspace, 

W  =  K  UD+ 1  •”  U^-|J  (24) 

where  u^Up...^,  are  the  ordered  eigenvectors  of  Rnw  arranged  from  the  largest  corresponding 

eigenvalues  to  the  smallest.  The  product  WWw  is  the  projection  operator  that  projects  the  array  response 
vector  onto  the  noise  subspace. 
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5.2.3  Through-the-Wall  Radar  Data 

In  through-the-wall  radar,  the  antennas  are  not  typically  placed  against  the  wall,  hut  rather  at  a  standoff 
distance.  This  is  mainly  due  safety  and  logistic  operation  reasons.  In  this  case,  and  for  a  single  wall,  the 
signal  travels,  on  transmission,  in  the  air,  through  the  wall,  and  then  in  the  air,  again.  It  propagates  with 
the  same  order  when  traveling  from  the  target  back  to  the  receiver.  The  signal  propagation  in  the  wall  is 
different  from  that  in  the  air  [16].  Consider  a  homogeneous  wall  with  known  dielectric  constant  £ .  The 
signal  path  from  the  antenna  to  a  position  through  the  wall  is  shown  in  Fig.  2.  Instead  of  following  a 
straight-line  direct  path  from  the  antenna  to  the  target  when  no  wall  is  present,  the  signal  now  propagates 
along  g, ,  g2 ,  and  g3.  From  Snell's  law,  the  angle  (pt  is 


<Pn  ~  sin  1 


sin  6n  1 

VF 


(25) 


If  the  targets  are  behind  a  wall,  the  angle  0n  as  well  as  the  range  to  the  center  rn  in  equation  (7)  should  be 
replaced  by  0n  and  rn  such  as. 


rn  =Sl+g2+Sv 


(26) 


Fig.  2  Geometry  of  through-the-wall  radar.  The  dotted  line  represents  direct  path  from  the  antenna  to  the 
center  of  the  scene.  However,  signal  travels  through  the  solid  line. 

The  parameters  g. ,  i  =  1,2,3  and  0n  can  be  calculated  from  the  following  equations  [17] 
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(27) 


8 .  =  -^a>/cos^> 

(xc  -  (xn]  -  yla)  tan  K  ))2  +  y2c  =  gl+gl-2gTig2co*(K+(pn-en\ 

where  (jc^0),^fl>)  is  the  position  of  the  /2-th  antenna.  The  above  equations  can  be  solved  numerically  by 
the  Newton  method  with  reasonable  initial  values.  Once  r  is  obtained  for  all  /z,  the  phase  term  can  be 
compensated,  resulting  in  equation  (7)  with  6n  in  place  of  6n.  Extension  of  equation  (27)  can  be  devel¬ 
oped  in  the  case  where  the  target  is  behind  multiple  walls. 

5.3  Beamspace  MUSIC  for  Imaging 

5.3.1  Beamspace 

As  stated  in  the  previous  section,  the  MUSIC  algorithm  can  be  applied  upon  generating  the  2-D  uniformly 
sampled  spatial  spectrum.  An  alternative  to  data  interpolation  can  be  provided  though  the  use  of  beam¬ 
forming.  Consider  the  2-D  delay-and-sum  (DS)  spatial  beams, 

1  n(-\n-\ 

b(x,  y)  =  -----  X  X  vv(m,/Oz(m,/i)exp  {j4nfmR„/c}  (28) 

A  m=r0  /i=() 

where  Rn  is  the  range  between  the  n- th  antenna  and  the  beamforming  point  (jc,y),  and  w(m,n)  is  an 

optional  window  that  can  be  used  to  shape  the  point  spread  function  in  terms  of  mainlobe  and  sidelobe 
characteristics.  If  a  wall  is  present  between  targets  and  antennas,  the  range  Rn  should  be  calculated 
through  equations  (27-28).  The  beams,  corresponding  to  uniformly  spaced  points  along  the  x-axis  and  y- 
axis,  produce  a  2-D  image  that  is  same  as  the  one  produced  by  the  backprojection  method  [18].  The 
M  xN  equally  spaced  beams  are  constructed  for  xt  =  X0  +  /A*  and  y*  =  Y()  +£Ay  ,  where  Av  and  Ay 
are  the  distances  between  beams  along  x-axis  and  y-axis,  respectively.  We  can  express  the  beams  using 
simpler  notations  as, 

b{l,k):=b(x,,yk)  (29) 
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for  /  =  0,...,Af  -1  and  k  =  0,...,/V-l  .  The  two-dimensional  discrete  Fourier  transform  (DFT)  of  b(l,k) 


is, 


B(m,n )  =  — —  YdYP^k^exP^~J(mlkx  +  nkky)) 


(30) 


=  ^ZWp>kf)e*PHW,kx+nkpky)) 


for  m  =  -1  and  n  =  0,...,  N  -1  .  In  the  above  equation,  kx-2n!M,  ky  =  2n  /  /V  ,  and  ( lp,kp )  is 

the  pixel  index  of  the  p  -th  target.  Note  that  b{l,k)  =  0  when  there  is  no  target  at  (/,&) .  Since  (30)  has  the 
same  structure  as  (11),  we  can  proceed  with  the  MUSIC  estimator,  without  the  need  for  performing 
interpolations.  In  this  case,,  the  sub-matrix  F.  .  is  redefined  as 


B(i,j  +  L- 1)  1 

B{i  + 1 ,  y  +  L  —  1 ) 


B(iJ)  S(i,y  + 1) 

B(i  + 1 ,  y)  +  1,7  +  1) 


5(/  +  AT  —  1 , 7)  ^( i  K  — 1,7  +  1)  *  *  *  B(i  +  K  —  1 , 7*  +  L*  ~  1 )  J 


and  the  covariance  matrix  follows  the  same  process  described  in  Section  5.2.2. 

5.3.2  Advantages 

Although  the  BS-MUSIC  requires  additional  computations  due  to  DS  beamforming,  it  has  key  advantages 
over  the  ES-MUSIC.  First,  it  is  robust  to  low  SNR.  The  ES-MUSIC  uses  2-D  spectrum  which  is  an 


interpolated  data  set  z(m,n)  of  the  received  signal.  Since  it  employs  raw  data,  the  target  returns  will  be 


highly  corrupted  by  noise  at  low  SNR,  resulting  in  a  compromised  quality  image.  In  contrast,  the  BS- 
MUSIC  enjoys  the  processing  gain  of  the  beamforming.  It  is  noted  that  the  beamforming  gain  is  up  to 
NsNf  when  a  target's  location  matches  the  beamforming  point  [19].  Second,  the  BS-MUSIC  does  not  put 

any  demands  on  target  characteristics.  The  ES-MUSIC  high-resolution  image  formation  process  begins 
with  a  very  important  assumption  of  point  targets.  Under  this  assumption,  the  signal  model  of  (1)  be¬ 
comes  valid.  However,  in  practice,  both  extended  targets  and  point  targets  could  be  present  in  the  same 
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scene.  For  instance,  interior  walls,  closets,  bookshelves,  large  and  small  file  cabinets,  and  tables  of 
different  sizes  are  typically  encountered  in  through-the-wall  radar  imaging  [20]. 

Consider  an  extended  target  located  at  (xe,ye) .  The  width  and  length  of  the  target  are  denoted  as  dx 
and  dv ,  respectively.  Assume  that  the  reflection  coefficient  cr  is  constant  over  the  target  2-D  extent. 
Then,  the  received  signal  becomes 

dj 2  r  4  r  } 

zext{myn)  -  a  f  exp<{ — —cos 0H(xe+x)\dx 

J dxn  [  c  J 


dj  2 


•  j  exp!j 

J~dyl2 


I  4^/W 


sin<9,(^+y)r^ 


=  sine 


2 7Ufm  COS  6dr 
J  m  n  x 


sine 


2nfm  sin  0ndx 


(31) 


cr4dxdv  exp| cos 6n  -  yr  sin  <9„ ) 


where  sinc(-)  is  the  Sine  function  given  by, 


sine(x)  = 


sin  x 


The  two  multiplicative  Sine  terms  in  (31),  due  to  non-zero  dx  and  dy ,  make  the  noise-free  virtual 
snapshots  zj ;  .  independent  of  each  other.  The  uniformly  sampled  version  of  zext{myn)  takes  the  form, 

'2 nkndy  ' 


„  ,  ,  .  2nk  d 

Zext\m'n)  =  SinC 


■cr4dxdv 


v  c 

(.  -4n 


sine 


expj -j—{kmX'-knyr) 


(32) 


Let 


u(m,n)  =  sine 


Ink  d  \  . 

2nknds  \ 

- 2L-JL  1  s  me 

C  J 

l  c  J 

(33) 


Then, 


("'*«)  =  <r4dxdyu(m,n)exp^-j^-(kmxe 


(34) 
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Note  that  zext(m,n)  =  z(m,n)  when  dx  =  dy  =  0 .  In  order  for  the  MUSIC  estimator  to  work,  the  noise-free 

virtual  snapshots  should  lie  in  the  signal  subspace  and  be  similar  in  structure  to  the  array  manifold  vector 
a(jt,y)  in  equation  (23).  However,  in  the  underlying  case, 

z.j=\J(i,j)a(xe,yt),  (35) 

where  U (i9j)  is  a  KLxKL  diagonal  matrix  such  that, 

diag[V(i,  j )}  =  [u{i,  j)  u(i  +  I  u{i  +  K  - 1 ,  j  +  L-\  )f  (36) 

The  values  of  u(i,j)  are  not  given,  as  they  are  determined  by  the  target  unknown  spatial  extent.  Further¬ 
more,  the  virtual  snapshots  for  the  same  target  are  different  since  the  elements  of  U (i9j)  vary  depending 
on  the  index  i  and  j .  This  difference  is  more  than  just  a  phase  change.  This  means  that  the  virtual 
snapshots  are  all  independent  (or  dissimilar),  and  subsequently,  the  rank  of  the  noise-free  covariance 
matrix  estimate  will  exceed  the  number  of  targets.  Accordingly,  applying  MUSIC  estimator  w  ith  the  array 
manifold  vector  a(jt,y)  fails  and  becomes  inappropriate  for  extended  targets. 

For  the  BS-MUSIC,  we  consider  the  beam  b(l,k)  to  be  within  the  extended  target.  Then, 

v  2  v  2 

b(l,k)  =  a  [  [  f(x\y)dxdy\  (37) 

-\/2  J  V2 

where  / (;t,  y)  represents  the  beam  pattern  of  the  beamformer,  and  bx  and  by  are  the  beamwidth  along  x- 

axis  and  y-axis,  respectively.  Suppose  that  an  extended  target  encompasses  the  beams 
b(l{)  —m{)9k0  - n0),b(l{) -m0  +  1,&{) -n0),...b(/0  +m0,k{)  +h0)  ,  as  shown  in  Fig.  3.  Note  that  the  center 

location  is  now  indexed  as  (/0,&0)  and  the  size  is  2m()x2n0.  The  corresponding  2-D  discrete  Fourier 
transform  of  the  beamspace  is, 

/0+mo  *0+no 

B(m,n)=  Yj  X  b(l>k)exp[-j(mkj  +  nky.k)}.  (38) 

If  b(l,k)  is  constant,  the  above  equation  becomes. 
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Fig.  3  Extended  targets  in  BS-MUSIC. 


B(m,n)  =  b(l0,k0)ex.p{-j(mkJQ -nkyk0)} 

"o 

■X^{mkx(lQ+l)}  (39) 

/=! 

n0 

•  2]  sm{nky(k0  +k)}, 

k  =  \ 

which  has  the  same  shortcoming  as  equation  (34).  However,  the  heampattern  of  an  array  f(xyy) 
changes,  depending  on  the  beamforming  point.  We  note  that  the  beamwidth  varies  as  a  function  of 
standoff  distance  and  the  imaging  point  relative  to  the  boresight  direction.  Therefore,  if  we  redefine 
b(l,k)  as. 


then. 


b(l,k)  =  <7  Jj/  (x,y\l,k)dxdy. 


b{l,k)*b{l\k\ 


for  k  ^  k'  or  /  ^  /' .  Therefore,  the  2-D  spectrum  of  beamspace  remains  as 


Z^+mo  *o+rto 

B(m,n)  =  z  z  b(lyk)exp{-j(mkj  +  nkyk)}. 

l=li)-mi)k=k0-n<) 


(40) 


(41) 
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The  above  equation  has  the  same  structure  as  equation  (30),  since  the  phase  terms  from  the  targets  do  not 
collapse  into  the  Sine  term.  The  extended  target,  in  essence,  is  viewed  as  an  aggregate  of  consecutive 
point  targets,  allowing  the  MUSIC  estimator  to  function  properly. 


5.3.3  Forming  Beamspace 
5.3.3. 1  Number  of  Beams 

It  is  assumed  that  all  antennas  of  the  imaging  system  have  identical  characteristics.  The  beamforming 
positions  can  be  located  anywhere  within  the  antenna  radiation  pattern.  Accordingly,  the  number  of 
beams  can  be  chosen  arbitrary.  However,  since  the  beams  must  be  uniformly  spaced  for  DFT  implemen¬ 
tations,  the  number  of  beams  is  decided  based  on  beam  spacing.  The  resolution  offered  by  the  imaging 
array  through  the  corresponding  DS  beamformer  presents  a  good  basis  to  form  the  different  beam 
orientations.  In  order  to  ensure  that  all  targets  in  the  scene  are  illuminated,  the  beams  should  be  separated 
by  less  than  a  beamwidth,  which  also  defines  target  resolution.  Target  resolution  includes  the  down-range 
resolution  and  the  cross-range  resolution.  The  down-range  is  the  range  perpendicular  to  the  array  plane, 
whereas  the  cross-range  is  parallel  to  the  array  plane.  The  down-range  resolution  of  the  DS  beamformer  is 
given  by 


where  B  is  the  signal  bandwidth.  The  cross-range  resolution  is 


~  _  rA 

a  ~2L 


(42) 


(43) 


where  r  is  the  down-range,  A  is  the  signal  wavelength,  and  L  is  the  array  length  [17].  If  the  beam 
spacing  is  chosen  as  the  beamwidth,  the  targets  positioned  close  to  or  at  the  middle  of  two  consecutive 
beams  will  be  illuminated  with  small  SNR.  Therefore,  it  would  be  better  if  the  beams  are  formed  more 
densely.  In  the  simulations,  we  used  half  the  cross-range  resolution  (one  fourth  of  the  beamwidth)  as  the 
beam  separation  (spacing). 


164 


5.3.3.2  Covariance  Matrix  Estimation 


When  the  distance  between  beams  are  less  than  the  beamwidth,  the  presence  of  a  point  target  will  be 
reflected  in  several  consecutive  beams,  leading  to  broadening  of  the  target  extent  and  representation  in  the 
beamspace.  In  this  case,  the  spatial  frequency  bandwidth  of  the  target  is  reduced  by  the  fixed  time- 
bandwidth  product.  Consequently,  the  2-D  spectrum  of  the  target  has  a  narrowband-like  representation. 
Therefore,  only  part  of  the  DFT  of  the  beams  should  be  included  in  covariance  matrix  estimation.  Let  Tx 
be  the  beamwidth  of  the  DS  beamformer  and  A tx  be  the  space  between  beams,  both  along  the  x-axis. 
Suppose  that  both  parameters  have  the  following  relation. 

A  tx=Tx/q.  (44) 

Then,  the  support  of  beams  (spatial  frequency  bandwidth)  will  be  1  /  (qAtx)  and  the  number  of  rows  of  2- 
D  spectrum  that  can  be  used  in  the  covariance  matrix  estimation  is 

M  ,ff  =  min  {A/,  La/  (45) 

w  here  L  J  represents  the  maximum  integer  that  does  not  exceed  the  number.  By  the  same  argument, 
when  A ty  =  rylp,  the  number  of  columns  that  include  spectrum  of  the  point  target  is 

Nfff  =  mm{N,lN/  p]},  (46) 

where  Tx  and  A tx  are  the  beamwidth  of  the  DS  beamformer  and  beam  spacing,  respectively,  both  along 
the  y-axis. 

When  q  and  p  are  larger  than  1,  implying  that  the  distance  between  beams  is  smaller  than  the 
beamwidth,  we  can  use  only  M eff  rows  and  Neff  columns  in  the  covariance  estimation  (see  Fig.  4).  The 

location  of  the  columns  and  rows  can  be  any  where,  depending  on  the  phase  of  the  reflection  coefficient  of 
the  target.  In  order  to  locate  the  spectrum  of  the  target,  we  used  the  magnitude  of  the  beams  prior  to  DFT 
implementation. 
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Fig.  4  Effective  elements  in  the  2-D  spatial  spectrum.  Among  M  xN  2-D  spectrum  of  the  beams,  only 

the  shaded  part  has  information  on  the  targets. 


b<ytn3»2ed6pa(5JS»^«ncy  Ncfro*ze<i  tfutai  tequtncy 

(a)  With  Phase  (b)  Magnitude  Only 

Fig.  5  Spatial  spectrum  of  the  DS  image  (a)  with  phase  information  and  (b)  with  magnitude  only.  The 
information  on  target  location  is  partly  present  and  it  is  located  around  the  center  when  only  magnitude  of 

the  image  is  used. 

In  some  applications,  such  as  target  recognition,  phase  information  from  target  may  be  necessary. 
However,  in  high-resolution  image  formations,  the  magnitudes  of  the  beams  are  sufficient.  After  exclud¬ 
ing  the  phase,  the  signal  is  located  at  the  center  of  the  2-D  spectrum,  and  the  center  part  in  the  covariance 
estimation  can  be  considered  (See  Fig.  5).  Extensive  simulations  have  shown  that,  depending  on  SNR, 
there  could  be  significant  difference  in  the  BS-MUSIC  performance  when  considering  the  entire  extent  of 
the  2-D  spectrum  vs.  only  using  the  region  where  most  of  the  signal  power  is  concentrated. 
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5.4  Simulations 


The  proposed  algorithm  is  tested  with  real  data  as  well  as  synthesized  data. 

5.4.1  Experiment 

Through-the-wall  data  set  is  collected  at  the  test  site  of  the  Center  for  Advanced  Communications  (CAC) 
at  Villanova  University,  PA,  USA  [20].  The  RF  system,  used  for  data  collection,  is  housed  in  the  Radar 
Imaging  Lab  of  the  CAC.  It  consists  of  off-the-shelf  equipment,  including  the  network  analyzer 
(ENA5071B,  Agilent,  Santa  Clara)  and  horn  antennas  (BAH-H1479,  BAE  Systems).  The  network 
analyzer,  which  has  operational  frequency  range  of  300kHz  -  8.5  GHz,  was  used  over  the  limited  fre¬ 
quency  band  2GHz  -  3GHz.  The  transmission  power  was  5dBm.  The  wideband  horn  antenna,  whose 
operation  frequency  range  is  1  -  12.4  GHz,  were  used  for  both  transmission  and  reception.  In  order  to 
move  the  antenna  and  position  it  precisely  over  the  P  different  locations,  the  Damaskos  Inc.  Field  Probe 
Scanner  model  7X7Y  was  used.  Although  the  signals  were  transmitted  and  received  on  the  2-D  vertical 
plane,  only  data  collected  at  the  horizontal  linear  position  is  used  for  the  simulation  Fig.  6  shows  the 
block  diagram  of  the  data  collection  system.  A  12'  x  8’  typical  interior  residential  wall  stood  between  the 
antenna  and  the  scene.  The  wall  consists  of  exterior  grade  3/4  inch  plywood  on  one  side  and  5/8 -inch 
gypsum  wallboard  was  applied  to  the  other  Between  them,  2x4  wood  studs  placed. 

Two  scenes  were  used  for  data  collection.  In  the  first  scene,  trihedrals,  dihedrals,  cylinder  and  sphere 
are  positioned  as  point  targets  in  a  25’  x  25’  room  (see  Fig.  7).  Reflections  from  interiors  other  than  targets 
were  minimized  by  absorbing  wall  modules  utilizing  18”  pyramidal  foam  absorber  and  laminated  polyu¬ 
rethane  foam  sheet  absorber.  Step-frequency  narrowband  signals  between  2  GHz  and  3  GHz  with  5  MHz 
interval  are  transmitted  and  received  by  the  horn  antenna  at  each  location,  and  this  is  repeated  at  57 
uniformly  spaced  locations.  The  total  array  length  is  1.2446  m.  In  the  other  scene,  a  6-ft  (1.82  m)long  and 
2.25  inch  (5.7  cm)  thick  wall  is  used  as  an  extended  target.  In  this  data  set,  the  number  of  frequencies  is 
335  and  the  number  of  antenna  locations  is  8 1 .  The  array  length  is  1 .778  m. 
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In  order  to  compare  the  performance,  we  processed  synthesized  data  of  similar  settings  where  we  can 
control  SNR.  Three  point  targets  are  considered  in  the  synthesized  data.  Fig.  8  summarizes  the  specifica¬ 
tion  of  data  sets. 

5.4.2  Results 

In  all  simulations,  DS  beamformer  with  a  rectangular  window  is  used.  In  DS  beamforming,  the  space 
between  beams  is  one  fourth  the  down-range  resolution  {q,  p  -  4).  The  center  half  of  the  2-D  spectrum  is 
used  in  the  covariance  matrix  estimation.  The  number  of  pixels  in  the  final  image  rendered  by  the  MUSIC 
estimator  is  same  as  the  number  of  beams. 

Fig.  9-12  are  the  log-scaled  images  obtained  from  the  synthesized  data.  The  true  target  locations  are 
indicated  by  Three  well-separated  targets  were  considered.  The  targets  locations  in  cross-range  and 
down-range  were  (-0.9,  4.1),  (-0.5  4.5)  and  (0.3,  3.4)  meter.  The  size  of  the  2-D  spectrum  and  sub¬ 
matrices  is  set  to  160x160  and  5x5,  respectively.  The  SNR  for  each  antenna  is  -20  dB.  The  DS  beam¬ 
forming  in  Fig.  9  shows  the  three  targets.  However,  it  also  shows  many  target-like  objects  due  to  high 
sidelobes 
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Data  Acquisition 
Positioner  Controller 
Computer 


Fig.  6  Block  diagram  of  the  data  collection  system. 
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Fig.  7  Ground  truth  of  the  targets  in  the  point-target  experiment  all  in  meters.  Note  that  the  targets  are  not 

in  the  same  x-y  plane  due  to  height  differences. 


Criterion 

Experiment  1 

Experiment  2 

Synthesized 

Bandwidth 

1  GHz 

1  GHz 

1  GHz 

Array  length 

1 .2446  m 

1.778  m 

1.2  m 

Freq.  step 

5  MHz 

3  MHz 

5  MHz 

Antenna  space 

0.022  m 

0.022  m 

0.02  m 

Fig.  8  Data  specification. 


and  noise.  Because  of  low  SNR,  the  ES-MUSIC  estimator  fails  to  identify  the  three  point  targets,  whereas 
the  BS-MUSIC  estimator  in  Fig.  1 1  and  Fig.  12  succeed  in  resolving  them.  The  difference  in  performance 
between  the  ES-MUSIC  and  the  BS-MUSIC  underscores  the  advantage  of  the  SNR  gain  provided  by  DS 
beamforming.  Comparing  the  two  BS-MUSIC  images.  Fig.  12  shows  better  result  than  that  of  Fig.  11. 

This  is  expected  since  the  covariance  estimate  Rrt</vv  in  Fig.  12  is  superior  to  R  in  Fig.  1 1. 

Real  data  measurements  consisting  of  point  targets  are  also  processed.  In  this  case,  the  beamspace 
size  is  also  160x160,  but  the  size  of  the  sub-matrices  is  14x14  and  the  dimension  of  the  signal  subspace 
is  25.  In  the  real  data  case,  since  the  number  of  point  targets  is  larger  than  the  synthesized  data,  we  used 
higher  dimension  sub-matrices  and  signal  subspace.  Once  again,  the  output  of  BS-MUSIC  (Fig.  15) 
shows  a  better  image  compared  to  the  ES-MUSIC  (Fig.  14).  Note  that  there  is  bias  in  the  MUSIC  target 
location  estimates.  This  is  because  the  scene  is  a  horizontal  2-D  plane  whose  height  is  the  same  as  the 
antenna,  but  the  targets  are  located  at  different  heights.  The  height  difference  results  in  bias  in  the  image. 
In  Fig.  14,  only  three  targets  are  identifiable,  while  Fig.  15  shows  seven  out  of  nine  targets. 

Real  wall  data  was  collected  and  used  to  verify  the  advantage  of  the  BS-MUSIC  in  dealing  with  an 
extended  target.  This  time,  8x8  sub-matrices  are  used  and  the  dimension  of  the  signal  subspace  is  four. 
Fig.  16  is  the  output  of  DS  beamformer.  Figure  17  and  Fig.  18  show  the  results  of  the  ES-MUSIC  and  the 
BS-MUSIC,  respectively.  The  image  by  the  ES-MUSIC  looks  meaningless  while  the  BS-MUSIC  clearly 
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shows  wall  image.  The  experiment  and  imaging  results  support  the  aforementioned  discussions  which 


consider  the  ES-MUSIC  as  a  good  imaging  approach  for  extended  targets. 

5.4.3  Conclusions 

High-resolution  image  formation  is  important  for  ISAR  and  TWI  applications.  The  ES-MUSIC  is  com¬ 
monly  used  to  obtain  high-resolution  imaging.  We  propose  the  BS-MUSIC,  which  applies  the  MUSIC 
algorithm  to  the  2-D  spatial  spectrum  in  beamspace.  The  key  advantages  of  the  BS-MUSIC  over  the  ES- 
MUSIC  are  its  robustness  to  low  SNR  and  ability  to  resolve  extended  targets.  These  two  advantages  make 
the  BS-MUSIC  more  suitable  for  urban  sensing  applications,  especially,  TWI,  where  a  typical  scene  is  a 
mixture  of  both  extended  targets  and  point  targets  with  several  reflections  from  surrounding  scatterers  are 
present.  The  performance  of  the  BS-MUSIC  with  both  synthesized  and  real  data  was  shown  to  be  superior 
to  the  ES-MUSIC. 
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Cross  range,  meter 

Fig.  9  The  target  locations  are  marked  by  s+\  Although  three  targets  are  visible,  many  false  small  targets 

appear  due  to  low  SNR  (-20  dB). 
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Fig.  10  Poor  imaging  quality  is  due  to  the  low  SNR  at  each  antenna. 
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Fig.  11  BS-MUSIC  using  the  covariance  matrix  estimated  with  diagonal  sub-matrices  as  in  eq.  (20). 
Although  three  targets  are  identifiable,  it  shows  artifacts  and  the  point  targets  are  blurred. 
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Fig.  12  BS-MUSIC  using  the  covariance  matrix  estimate  in  eq.  (21).  Note  that  sharper  image  is  obtained 
compared  to  the  DS  image. 
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Fig.  13  DS  beamformer  for  the  experimental  data  (point  targets). 
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Fig.  14  ES-MUSIC  for  experimental  data.  Only  three  targets  out  of  nine  targets  are  identifiable. 


Fig.  15  BS-MUSIC  for  the  experimental  data.  Seven  out  of  nine  targets  are  identifiable.  Note  that  the  bias 
is  due  to  the  fact  that  the  heights  of  imaged  targets  are  different  from  the  2-D  plane. 
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Fig.  16  DS  beamformer  for  the  experimental  wall  data. 
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Fig.  17  BS-MUS1C  for  the  experimental  wall  data.  The  BS-MUSIC  fails  to  image  the  wall. 
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Fig.  18  BS-MUSIC  for  experimental  wall  data.  A  wall  is  clearly  seen  in  this  image. 
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